{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1D viscoelastic SH modelling\n",
    "\n",
    "After deriving the equations of motion for 1D wave propagation in viscoelastic SH media, we can now solve the problem using a staggered grid FD scheme. Furthermore, the anelastic coefficients of the Generalized Maxwell body (GMB) are optimized by a global Differential Evolution (DE) algorithm to achive a constant Q-spectrum. Finally, elastic and viscoelastic modelling results for a homogeneous model are compared."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FD solution of 1D isotropic viscoelastic SH problem\n",
    "\n",
    "As derived in the last lesson, we can describe the 1D viscoelastic SH problem by the partial differential equations:\n",
    "\n",
    "\n",
    "\n",
    "The 2nd order space/time FD scheme on a staggered grid \n",
    "\n",
    "<img src=\"images/SG_1D_SH-Cart.png\" width=\"50%\">\n",
    "\n",
    "using explicit time-stepping with the Leapfrog method to solve the 1D viscoelastic SH problem can be written as \n",
    "\n",
    "\\begin{align}\n",
    "\\scriptsize \n",
    "\\frac{v_y(i,n+1/2) - v_y(i,n-1/2)}{dt} &\\scriptsize= \\frac{1}{\\rho}\\biggl\\{\\biggl(\\frac{\\partial \\sigma_{yx}}{\\partial x}\\biggr)^c(i,n) + f_y(i,n)\\biggr\\}, \\notag\\\\\n",
    "\\scriptsize\\frac{\\sigma_{yx}(i+1/2,n+1) - \\sigma_{yx}(i+1/2,n)}{dt} &\\scriptsize= \\mu_{x,u} \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,n+1/2) - \\sum_{l=1}^L Y_l\\xi_l(i+1/2,n),\\notag\\\\\n",
    "\\scriptsize \\frac{\\xi_l(i+1/2,n+1/2) - \\xi_l(i+1/2,n-1/2)}{dt} &\\scriptsize= \\omega_l \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,n+1/2) - \\omega_l \\xi_l(i+1/2,n-1/2)\\notag\n",
    "\\end{align}\n",
    "\n",
    "with the unrelaxed shear modulus \n",
    "\n",
    "\\begin{equation}\n",
    "\\mu_u = \\frac{\\mu}{1-\\sum_{l=1}^L Y_l}\\notag\n",
    "\\end{equation}\n",
    "\n",
    "and the harmonically averaged unrelaxed shear modulus\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu_{x,u} = 2 \\biggl(\\frac{1}{\\mu_u(i)}+\\frac{1}{\\mu_u(i+1)}\\biggr)^{-1}\\notag\n",
    "\\end{equation}\n",
    "\n",
    "Notice that we have to evaluate the memory variables $\\xi_l$ in the stress update at a full time step, which has to be estimated from the $\\xi_l$-values at half-time steps by arithmetic averaging:\n",
    "\n",
    "\\begin{equation}\n",
    "\\xi_l(i+1/2,n) = \\frac{\\xi_l(i+1/2,n+1/2)+\\xi_l(i+1/2,n-1/2)}{2}\\notag\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial and boundary conditions\n",
    "\n",
    "Because we have analytical solutions for wave propagation in homogeneous elastic media, we should test our code implementation for a similar medium, by setting density $\\rho$ and shear modulus $\\mu$ to constant values $\\rho_0,\\; \\mu_0$\n",
    "\n",
    "\\begin{align}\n",
    "\\rho(i) &= \\rho_0 \\notag \\\\\n",
    "\\mu(i) &= \\mu_0 = \\rho_0 V_{s0}^2\\notag\n",
    "\\end{align}\n",
    "\n",
    "at each spatial grid point $i = 0, 1, 2, ..., nx$ and staggered grid points $i+1/2$ in order to compare the numerical with the analytical solution. For a complete description of the problem we also have to define initial and boundary conditions. The **initial condition** is \n",
    "\n",
    "\\begin{equation}\n",
    "v_y(i,-1/2) = \\sigma_{yx}(i+1/2,0) = \\xi_{l}(i+1/2,-1/2) = 0, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "so the modelling starts with zero particle velocity, shear stress and memory variable amplitudes at each spatial grid point. As **boundary conditions**, we assume \n",
    "\n",
    "\\begin{align}\n",
    "v_y(0,n) &= \\sigma_{yx}(1/2,n) = \\xi_{l}(1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(nx,n) &= \\sigma_{yx}(nx+1/2,n) = \\xi_{l}(nx+1/2,n) = 0, \\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "for all full and staggered time steps n, n+1/2. This **Dirichlet boundary condition**, leads to artifical boundary reflections which would obviously not describe a homogeneous medium. For now, we simply extend the model, so that boundary reflections are not recorded at the receiver positions.\n",
    "\n",
    "Let's implement it ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "from numba import jit\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "from scipy.optimize import differential_evolution\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimizing Yl coefficients\n",
    "\n",
    "To achieve a constant Q-spectrum, we have to optimize the anelastic coefficients $Y_l$. This can be achieved by minimizing the objective function:\n",
    "\n",
    "\\begin{equation}\n",
    "\\chi(Y_l) = \\int_{\\omega_{min}}^{\\omega_{max}} \\biggl(Q^{-1}(Y_l,\\omega) - Q^{-1}_{opt}\\biggr)^2 d\\omega\\notag\n",
    "\\end{equation}\n",
    "\n",
    "where $Q^{-1}_{opt}$ denotes the target constant inverse Q-value within the frequency range of the source wavelet and \n",
    "\n",
    "\\begin{equation}\n",
    "Q^{-1}(Y_l,\\omega) = \\sum_{l=1}^L Y_l \\frac{\\omega_l\\omega}{\\omega^2+\\omega_l^2}\\notag\n",
    "\\end{equation}\n",
    "\n",
    "is an approximation of the GMB Q-spectrum for $Q>>1$ according to ([Blanch et al. 1995](https://library.seg.org/doi/abs/10.1190/1.1443744), [Bohlen 2002](https://www.sciencedirect.com/science/article/pii/S0098300402000067), [Yang et al. 2016](https://jean-virieux.obs.ujf-grenoble.fr/IMG/pdf/2016_yang_viscoelastic_final-3.pdf)). The objective function can be minimized using local or global optimization methods. In this case, we will use a global optimization [Differential Evolution (DE)](https://pablormier.github.io/2017/09/05/a-tutorial-on-differential-evolution-with-python/) algorithm available from `SciPy Optimization` library.\n",
    "\n",
    "In the following example, we want to solve the viscoelastic SH-problem by the (2,2)-FD algorithm defined above for a homogeneous model using 4 Maxwell bodies and target for a constant Q=20 value in a frequency band between 5 and 100 Hz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define GMB models\n",
    "# -----------------\n",
    "Qopt = 20     # target Q-value\n",
    "L = 4         # number of Maxwell bodies    \n",
    "nfreq = 50    # number of frequencies to estimate Q-model\n",
    "\n",
    "fmin = 5.     # minimum frequency \n",
    "fmax = 100.0  # maximum frequency\n",
    "\n",
    "f = np.linspace(fmin,fmax,num=nfreq)  # calculate frequencies to estimate Q-model\n",
    "w = 2 * np.pi * f                     # circular frequencies to estimate Q-model \n",
    "\n",
    "fl = np.linspace(fmin,fmax,num=L)  # calculate relaxation frequencies\n",
    "wl = 2 * np.pi * fl                # circular relaxation frequencies "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define the objective function to optimize the $Y_l$ values ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Objective function to optimize Yl values\n",
    "# ----------------------------------------\n",
    "def obj_Yl(Yl):                   \n",
    "    \n",
    "    # Calculate Qs model based on GMB\n",
    "    # -------------------------------\n",
    "    Qinv_GMB = np.zeros(nfreq)\n",
    "    Qinv_const = (1/Qopt) * np.ones(nfreq)\n",
    "    \n",
    "    for l in range (0,L):    \n",
    "        Qinv_GMB += Yl[l] * (wl[l] * w) / (w**2+wl[l]**2)\n",
    "    \n",
    "    # Calculate objective function\n",
    "    obj_Qinv = np.sum((Qinv_GMB - Qinv_const)**2)\n",
    "    \n",
    "    return obj_Qinv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A function to evaluate the Q-spectrum for given $Y_l$ values might be also useful ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate Q-spectrum for given Yl parameters and circular relaxation frequencies\n",
    "# --------------------------------------------------------------------------------\n",
    "def calc_Q(Yl):            \n",
    "    \n",
    "    # Calculate Qs model based on GMB\n",
    "    # -------------------------------\n",
    "    Qinv_GMB = np.zeros(nfreq)    \n",
    "    for l in range (0,L):    \n",
    "        Qinv_GMB += Yl[l] * (wl[l] * w) / (w**2+wl[l]**2) \n",
    "    \n",
    "    return f, Qinv_GMB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To optimize the $Y_l$ values, we define their bounds necessary for the DE-algorithm, minimize the objective function, store the resulting $Y_l$ values and plot the optimized Q-spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Final obj_Qinv =  7.651346890187525e-05\n",
      "Yl =  [0.08706609 0.03667265 0.         0.0629249 ]\n",
      "wl =  [ 31.41592654 230.38346126 429.35099599 628.31853072]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAFOCAYAAAD5H3jwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmUHXWd9/H3J93ZOwuQ0CQhEHhERXEM0HBQxrEDyIOOIzrHjccNRSPqOO7KOGeewWVcHhm3UWdEUXBGDQzLiIijCLSIsphAZMAoCkRIAlkgIXQ20sn3+eNXl77d6aXS6a7b6d/ndU6dW1W37q3v/XX1/dxfVd26igjMzMxyMa7RBZiZmVXJwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVhx8ZmaWFQefDYmkT0t67xAed7ukZw9w/zMk3SnpCUl/u29V7n8krZR0Wt30PZLah3kdF0v65HA+Z9X6aKce0xXWUbotG1Wj7cnBZ3uQNFXSJyXdVwTQbyW9ve7+2cAbga8P4ekvAD4+wP0fBjoiYlpEfHkIz78HSa+VdJukLZLWFePvlKTheP6RFBHPjoiOqtZXvDk/KWlWr/nLJYWkBVXVYjZSHHzWg6QDgJuBI4BTgenA24BPSDqnWOxs4NqI2DaEVVwNLJI0p5/7DwfuGcLzIqm5j3kfAL4EfA44BGgFzgVOBiYMZT1D1Vd9o9QDwFm1CUnPASY3rhyz4eXgs96+BKwFXh8RKyP5JfB54F3FMi8Gfl7/IEnTJK2W9MJe8+cXPYWDACJiO7AMOL33iiXdACwCviKpU9LTJR0tqUPSpmK338t6PWalpI9IugvYUh8ukmaQepfvjIjLI+KJ4vXcGRGvi4gddcvOlXSFpPWSHqjfzVqs44OS7pL0uKRLJU3ai8f2qE/Seb1606/o749R2z0m6TVFm9SGHZI6Blt/cf+xku4o1ncpMKmvddX5d1KPvuZNwHf6qK3P1yHpf0l6TNJxdfVtkNQu6c2Sflj3HH+UdFnd9EOSFpZ5XWUVbfih4u+3RdJFklol/bio/WfFBz5KbG/9tuVw1WsViAgPHogISL2tXcBxfdz3KuCxYnw9cEKv+88HftLP83YC7XXTXwY+38+yHcBbi/HxwB+Bj5J6Z6cATwDPqFt+JbAcmA9M7vVcZwBdQPMgr3scKYz/b7GeI4H7gf9dt47bgbnAgcAK4Ny9eGyP+oq2nFs89jXAFmBO3fKn9Xp9p/Wqd3pRw9tLrH8C8CfgfUV7vhLYCXyyn7ZYCZwG/B44GmgCHiq2jQAW9Nom+nsdbytqnAL8BLigmH8ksKl4zJyittV1920s7ivTrgO2U6/7biX19ucB64A7gGOBicANwD8yyPY2UFsOVu9gNXqodnCPz+qdBjwUEXf0cd88YFUxPpP0hgCApCbgHcA3i+nZko6se2wXPXeVPVE8x2BOAlqAz0TEkxFxA3ANdbvhCl+OiIdiz12vs4ANEdFVV+uvik/z2yT9RTH7BGB2RHy8WM/9wDeA1/Zax5qIeAz4IbBwLx/7VH0R8Z/Fc+2OiEuBPwAnlmgPJI0Dvkc6Dvr1Eus/ifQm/cWI2BkRlwO/LrGqWq/vRcDvgNW9FxjodUTEN4rp20gB9/fF/PtJf/+FwAtJobha0jOL6V9ExO4Sr2tv/UtErI2I1cAvgNsi9fx3AFeRQnCw7W2gthzuem0E7S/HHKwas+kOt95eQfpkDOlT+bS6+44BDia9iQG8HxBwnqTJxbLr6pafRvrUP5i5pCDeXTfvT6QQrvdQP49/FJglqbkWfhHxfABJq+je1X84MFdSfU1NpDfImkfqxrcWtZV9bI/6JL2R1EYLilktpJAu459I7VfbjTbY+ueSelT1vz/2pxLr+XfgJtKx3j12c0Kp1/EN0jHdxVG3W5m0m7wdeFoxvokUes+jexd6mXbdG2vrxrf1Md3C4NvbQG053PXaCHLwWb0HgMMljav/55f0IqANeH0x6y7g6XR/2p0HbIyIzcX0GaQ3TkhvaBuBO+vWczTwHyXqWQPM71XPYcC9vZbr70clbwF2AGcCVwywnoeAByLiqBI1DeWxT9Un6XBSIJwK3BIRuyQtJ31QGJCk15J6HydExM6S638YmCdJdW/YhwH3DbSuiPiTpAeAlwDn9L5/sNchqQX4InARcL6kK4reMqRw+ytSqH6KFHyvIwXfV0q+rpEw2PY2UFs2ol4bIu/qtHo/Km4/KWmKpImSXg98H3hVRNR6LteSAq3mMWC6pCMknUU6xvEsSTNJx/6+WHsjkTQROB64rkQ9t5GOG31Y0nil77P9FbCkzIuJiE3Ax4CvSXqlpBZJ44qTJ6bWLXo7sLk4CWWypCZJx0g6ocRq9vaxU0lBuB5A0ptJPeYBSToW+Bfg5RGxfi/WfwtpV/PfKp1Y89eU3K1KCrxTImLLEF7Hl4BlEfFW0nb1b3X3/Zx0EtPkiFhF6hWdARxE9wekffmbDNVg29tAbdmIem2IHHz2lIjoJH2Cfw7pQPw20q6sF0bEtXWLfgd4SbEbE1LPbwnpJI5zgJcBz6f7GM9n6x77MtLxqTUl6nmyWP7FwAbga8AbI+J3e/Ga/l/xGj5M2t26lvT9w48AvyqW2UV6g1tI6vVuIB2vnFHi+ffqsRHxW+CfSW+ia0lt/csSL+VM4ADgZnWf2fnjwdZftOFfk76CspF0EsqVJdZHRNwXEUv39nVIOpMUZOcWi78fOE7S64rH3ks64ekXxfRm0okgvyxezz79TYZqsO1toLZsRL02dOq5u9qsm6RXkz65P7tuN1Xtvk8B6yLii3v5nLcB50TE3cNXqZlZeQ4+G5CkdwD3RMRNja7FzGw4VB58xanvS0lnR71U0hGk3WQHkr5b84Zil4KZmdmwa8QxvveQvtha81ngC8XZUBvp4wwyMzOz4VJp8Ek6FPhLur/oLNLVES4vFrkEeHmVNZmZWV6q/h7fF0ln19W+/HwQsKnuyhqr2PPLyQBIWgwsBpg8efLx8+fPH+FS93+7d+9m3DifuDsYt1N5bqvy3FblDUdb3XvvvRsiYnaZZSsLPkkvJZ0FuEzdvy/W15d2+zzoGBEXAhcCtLW1xdKlfZ5lbXU6Ojpob29vdBmjntupPLdVeW6r8oajrSSVuSIRUG2P72TgZZJeQrqi+XRSD3Bm3SWlDiVdPcHMzGxEVNYPj4i/i4hDI2IB6cKtN0TE64AbSVc5h/TzJz+oqiYzM8vPaNgB/RHg/ZL+SDrmd1GD6zEzszGsIRepjogO0u+u1X6mpOy1A83MzPbJaOjxmZmZVcbBZ2ZmWXHwmZlZVhx8ZmaWFQefmZllxcFnZmZZcfCZmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVhx8ZmaWFQefmZllxcFnZmZZcfCZmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZaWy4JM0SdLtkn4j6R5JHyvmXyzpAUnLi2FhVTWZmVl+mitc1w7glIjolDQeuFnSj4v7PhQRl1dYi5mZZaqy4IuIADqLyfHFEFWt38zMDEApjypamdQELAOeBnw1Ij4i6WLgeaQe4fXAeRGxo4/HLgYWA7S2th6/ZMmSyureX3V2dtLS0tLoMkY9t1N5bqvy3FblDUdbLVq0aFlEtJVZttLge2ql0kzgKuDdwKPAI8AE4ELgvoj4+ECPb2tri6VLl454nfu7jo4O2tvbG13GqOd2Ks9tVZ7bqrzhaCtJpYOvIWd1RsQmoAM4IyIejmQH8G3gxEbUZGZmeajyrM7ZRU8PSZOB04DfSZpTzBPwcuDuqmoyM7P8VHlW5xzgkuI43zjgsoi4RtINkmYDApYD51ZYk5mZZabKszrvAo7tY/4pVdVgZmbmK7eYmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVhx8ZmaWFQefmZllxcFnZmZZcfCZmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVhx8ZmaWFQefmZllxcFnZmZZcfCZmVlWKgs+SZMk3S7pN5LukfSxYv4Rkm6T9AdJl0qaUFVNZmaWnyp7fDuAUyLiucBC4AxJJwGfBb4QEUcBG4FzKqzJzMwyU1nwRdJZTI4vhgBOAS4v5l8CvLyqmszMLD+KiOpWJjUBy4CnAV8FPgfcGhFPK+6fD/w4Io7p47GLgcUAra2txy9ZsqSyuvdXnZ2dtLS0NLqMUc/tVJ7bqjy3VXnD0VaLFi1aFhFtZZZt3qc17aWI2AUslDQTuAo4uq/F+nnshcCFAG1tbdHe3j5SZY4ZHR0duJ0G53Yqz21VntuqvKrbqiFndUbEJqADOAmYKakWwIcCaxpRk5mZ5aHKszpnFz09JE0GTgNWADcCrywWexPwg6pqMjOz/FS5q3MOcElxnG8ccFlEXCPpt8ASSZ8E7gQuqrAmMzPLTGXBFxF3Acf2Mf9+4MSq6jAzs7z5yi1mZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVhx8ZmaWFQefmZllxcFnZmZZcfCZmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVhx8ZmaWFQefmZllxcFnZmZZcfCZmVlWHHxmZpaVyoJP0nxJN0paIekeSe8p5p8vabWk5cXwkqpqMjOz/DRXuK4u4AMRcYekacAySdcV930hIi6osBYzM8tUZcEXEQ8DDxfjT0haAcyrav1mZmYAiojqVyotAG4CjgHeD5wNbAaWknqFG/t4zGJgMUBra+vxS5Ysqaja/VdnZyctLS2NLmPUczuV57Yqz21V3nC01aJFi5ZFRFuZZSsPPkktwM+Bf4qIKyW1AhuAAD4BzImItwz0HG1tbbF06dKRL3Y/19HRQXt7e6PLGPXcTuW5rcpzW5U3HG0lqXTwVXpWp6TxwBXAdyPiSoCIWBsRuyJiN/AN4MQqazIzs7xUeVangIuAFRHx+br5c+oWewVwd1U1mZlZfqo8q/Nk4A3A/0haXsz7KHCWpIWkXZ0rgbdXWJOZmWWmyrM6bwbUx13XVlWDmZmZr9xiZmZZcfCZmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVhx8ZmaWFQefmZllxcFnZmZZcfCZmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlpXmRhdgZmaZiYBNm2D1ali9mtaODmhvr2z1Dr7RbOtWWLUKHnoobSQ7d/Y9AEyZAi0tMHXqU8PU++6Dww6DWbNg2jSQGvt6zGzs274dHn44hdqaNWkoAu6pYc0a2LbtqYccDfAP/5DexyowaPBJOrDE8+yOiE3DUE9+duyApUvh1lvh/vtTyNWGRx/dp6c+oX5i4kQ4+OA0zJ7dPT53LsybB4cemoY5c2D8+H1ar5mNQbVAqw1r1ux5u2YNPPbY0J5/9Wo46qjhrbkfZXp8a4phoO5CE3DYsFQ01j3+OPzqV3DzzfCLX8Dtt6fwG2k7dnQH6kAkOOSQFIbz56ce4+GH97ydPdu9R7OxYPfuFFSPPJKGtWtTiD3ySHfA1cY3DXPfpqUlvc/Mm8cjTU0c0lzdDsgya1oREccOtICkO4epnrFp5Ur4znfgqqvgN79J+7fLaG5OvbD58+Ggg1JPbMKEdFs/QNot2tkJW7Y8NXSuXUvL7t2wbl2P3QoDiuje4Jcu7XuZSZO6g7CvYd68VLuZVW/XLtiwIYXY2rXp/7/+tjY88kia19U1vOtvakp7jmp7k+bO7R6vH6ZPf+ohv+vo4JAjjhjeOgZQ5t3ppGFaJi9bt8IVV8DFF8MNNwy87FFHwZ//OTznOd29rPnzobUVxg39xNulHR201w4Yb9kC69enDb02rF2bdk2sWpWG1avTP8Ngwbx9O9x7bxr60tSUNuzDDuvZU6wN8+fDjBlDfl1mWdm1K/XK1q/fc1i3bs//60cfLf/hem80N6f3pLlzu4Ntzpye4/PmpXMKmpqGf/3DqEzw3QIct6/LSJoPfAc4BNgNXBgRXyqOIV4KLABWAq+OiI0l6hp9ItJuzG9/Gy67DJ54Ys9lxo2DhQtT0L3gBen2kENGvrbaSS8LFgy83M6dqbf30EMpDP/0J3jwwZ63jz8+8HPs2pWWffDB/peZNq27N1s7vjh/fvcnxHnzUi/Xu1RtLNm1K+0yfPTRFGaPPpp6Z/VDbd769d3TIxFkNTNnpveg1tY01MLskEN6jh900D59EB9NygTf0ZLuGuB+AWU+vncBH4iIOyRNA5ZJug44G7g+Ij4j6TzgPOAjJZ5v9OjqSkH3qU/BPffsef+4cXD66XD22fDiF/fo4o8648d398z68/jj3SFYG1au7B5fu3bw9TzxBKxYkYb+TJjQc5dJ7R+wNrS2ptuDD/YJOVat7dth48YUYhs3puGxx9JQjB+9YgV87nM9Q27jxpENsZoDDugOstpw8MHd47X/oYMPTocuMlMm+J5ZYpldgy0QEQ8DDxfjT0haAcwDzgTai8UuATrYX4LvySfhP/4DPv1p+OMf97z/Gc9IYfeGN6Q37rFixoy0W/Y5z+n7/m3bUm+x1uur7zE++GDqTW7fPvh6nnyyO0wHc9BB6aSb/oZZs+DAA9Nw0EH+ekfOurrSB6/Nm9OHuM2bu8cffzyFWf1tbbw+6Epsv63DWfPMmX1v07VAqz9be9YsfxAchKKKTx+9VyotAG4CjgEejIiZdfdtjIgD+njMYmAxQGtr6/FLliypptg+jHvySQ750Y847NJLmdSrd9M1eTLrTj2VR844g83PelZD31w7OztpaWlp2Pr7FUHz5s1MXL+eSevXM7E2rFvHxA0bmPDoo0zcsIHmLVtGrITdTU10TZ/OzmnT2DFlCjFjBl1Tp9LV0tJzmDqVXVOmsGvKFLqK29p4TJgwYvWNVlVvU+rqomn7dsZt397jtmn7dpq2bUtD/fi2bTRt3ZqGbdtorhtv2ro1TZf50DWCuqZOZeeMGeycNi1tgzNm9D1Mn87OmTPZOX06McZPFhuO7WrRokXLIqKtzLKVB5+kFuDnwD9FxJWSNpUJvnptbW2xtL8zDkfSxo1w0UXwz/+cTgKpd8AB8J73wLvfnXoVo0BH/ckt+6MtW7q/H7R6dRqvnY1Wf/r1unXV7D7qbfz4HhcMYMqUPacnTYLJk/u+nTAhDRMn9hyvna3b3LznbXNzOnFg3Lg01MZrt7UPWn3dRvQcdu/uOb5rV//Dzp3Q1cWyW2/l+Oc+t/viCV1dqWfe37BjR+od1Ybe09u29T9s2TL8ZxwOl+bm9D9/wAGpN3bAAd17FIphxSOPcPTJJ3fvZajdN8ZDbCiG471KUungq/QvIGk8cAXw3Yi4spi9VtKciHhY0hxg3YgX8tnPpn+6006DE08cfLfAsmXwta/B97+/59cCZs+GD3wA3vGO0X3sbn80dSo87WlpGEhXVzp+Un+WW++hdoyldrt1677Xt3Nn2vU13N9vGsWOb3QBw0FK/6vTp6fd9rXx6dNTiM2YseftjBk9Q27KlEH35qzt6ODo/fmD5xhWKvgkzQVOBaYAv4uIn+/tiiQJuIj0vcDP1911NfAm4DPF7Q/29rn3SgR85SvpGNT556cvUb7whSkETzsNnv3stEFv355OWPnqV9OXzHubOxc+/GF429squ8yO9aN2mnXrXhxV2b79qRC888YbOfbII7tDrPcxntrxoCee6Dleu1ycjZympp696fpedUtL91Cbnjo1Hb8daGhpGTNnJ9rQlLlk2el0n3SyAzhX0hTgzRHxq71Y18nAG4D/kbS8mPdRUuBdJukc4EHgVXvxnHvv3ntT6NV0dsKPfpQGSG+eJ52UrqzS1yXDFi6Ed70rnbAyceKIlmojaNKkp75Y+/ijjw7tArk7dvS4YMAeQ21XXu/bbdvSY2u7Autva+NdXd1DbZdi7ba2W3L37p7ju4pzzGq7fetvI9IHutpQ2y1aP93U1P9Q7H7dvG0b02u762q7X3vvrq0fJk5MbT1pUvd4/bzJk/sfpk5Nz2E2zMr0+D4JvCAinjptUdLzgG8UYbUlIu4e7Eki4mb6v+zZqWWKHRZz5sD3vgc/+1kaen/XbO1a+EGvTueECfDqV8M735lC0WcDGqQ38IkTR80x3Srcsb8fNzajXPBNqA89gIi4RdJfA9eQeoH9nNc+Ck2fDmedlYYIuO++FIDXX5+usFJ/gdXDD0/H7t7ylnQsz8zM9ntlgm+7pNkRsb5+ZkTcK2kXVfbWhpvUffLEueem3UXLl6drVB5+eDrmN8ovvWNmZnunTPB9DvgvSa+KiDW1mZJmATsiYuTPwqzKuHFw3HFpMDOzMWnQ4IuIKyRNBG6RtAz4DTABeDXp+J+Zmdl+o9Q5vRHxPdKP5F5Dui7nTuD/RMQlI1ibmZnZsCv9BfaI2Ap8awRrMTMzG3H+FqeZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVhx8ZmaWFQefmZllxcFnZmZZcfCZmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVioLPknfkrRO0t11886XtFrS8mJ4SVX1mJlZnqrs8V0MnNHH/C9ExMJiuLbCeszMLEOVBV9E3AQ8VtX6zMzM+qKIqG5l0gLgmog4ppg+Hzgb2AwsBT4QERv7eexiYDFAa2vr8UuWLBn5gvdznZ2dtLS0NLqMUc/tVJ7bqjy3VXnD0VaLFi1aFhFtZZZtdPC1AhuAAD4BzImItwz2PG1tbbF06dIRrHRs6OjooL29vdFljHpup/LcVuW5rcobjraSVDr4GnpWZ0SsjYhdEbEb+AZwYiPrMTOzsa+hwSdpTt3kK4C7+1vWzMxsODRXtSJJ3wfagVmSVgH/CLRLWkja1bkSeHtV9ZiZWZ4qC76IOKuP2RdVtX4zMzPwlVvMzCwzDj4zM8uKg8/MzLLi4DMzs6w4+MzMLCsOPjMzy4qDz8zMsuLgMzOzrDj4zMwsKw4+MzPLioPPzMyy4uAzM7OsOPjMzCwrDj4zM8uKg8/MzLLi4DMzs6w4+MzMLCsOPjMzy4qDz8zMsuLgMzOzrDj4zMwsKw4+MzPLioPPzMyy4uAzM7OsOPjMzCwrlQWfpG9JWifp7rp5B0q6TtIfitsDqqrHzMzyVGWP72LgjF7zzgOuj4ijgOuLaTMzsxFTWfBFxE3AY71mnwlcUoxfAry8qnrMzCxPiojqViYtAK6JiGOK6U0RMbPu/o0R0efuTkmLgcUAra2txy9ZsmTkC97PdXZ20tLS0ugyRj23U3luq/LcVuUNR1stWrRoWUS0lVm2eZ/WVKGIuBC4EKCtrS3a29sbW9B+oKOjA7fT4NxO5bmtynNblVd1WzX6rM61kuYAFLfrGlyPmZmNcY0OvquBNxXjbwJ+0MBazMwsA1V+neH7wC3AMyStknQO8BngRZL+ALyomDYzMxsxlR3ji4iz+rnr1KpqMDMza/SuTjMzs0o5+MzMLCsOPjMzy4qDz8zMsuLgMzOzrDj4zMwsKw4+MzPLioPPzMyy4uAzM7OsOPjMzCwrDj4zM8uKg8/MzLLi4DMzs6w4+MzMLCsOPjMzy4qDz8zMsuLgMzOzrDj4zMwsKw4+MzPLioPPzMyy4uAzM7OsOPjMzCwrDj4zM8uKg8/MzLLi4DMzs6w0N7oAAEkrgSeAXUBXRLQ1tiIzMxurRkXwFRZFxIZGF2FmZmObd3WamVlWFBGNrgFJDwAbgQC+HhEX9rHMYmAxQGtr6/FLliyptsj9UGdnJy0tLY0uY9RzO5XntirPbVXecLTVokWLlpU9TDZagm9uRKyRdDBwHfDuiLipv+Xb2tpi6dKl1RW4n+ro6KC9vb3RZYx6bqfy3Fblua3KG462klQ6+EbFrs6IWFPcrgOuAk5sbEVmZjZWNTz4JE2VNK02DpwO3N3YqszMbKwaDWd1tgJXSYJUz/ci4r8bW5KZmY1VDQ++iLgfeG6j6zAzszw0fFenmZlZlRx8ZmaWFQefmZllxcFnZmZZcfCZmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWXHwmZlZVhx8ZmaWFQefmZllxcFnZmZZcfCZmVlWHHxmZpYVB5+ZmWXFwWdmZllx8JmZWVYcfGZmlhUHn5mZZcXBZ2ZmWRkVwSfpDEm/l/RHSec1uh4zMxu7Gh58kpqArwIvBp4FnCXpWY2tyszMxqqGBx9wIvDHiLg/Ip4ElgBnNrgmMzMbo0ZD8M0DHqqbXlXMMzMzG3bNjS4AUB/zYo+FpMXA4mKyU9LvR7SqsWEWsKHRRewH3E7lua3Kc1uVNxxtdXjZBUdD8K0C5tdNHwqs6b1QRFwIXFhVUWOBpKUR0dboOkY7t1N5bqvy3FblVd1Wo2FX56+BoyQdIWkC8Frg6gbXZGZmY1TDe3wR0SXpb4CfAE3AtyLingaXZWZmY1TDgw8gIq4Frm10HWOQdw2X43Yqz21VntuqvErbShF7nEdiZmY2Zo2GY3xmZmaVcfCZmVlWHHxjgKT5km6UtELSPZLeU8w/UNJ1kv5Q3B7Q6FpHC0lNku6UdE0xfYSk24q2urQ4wzh7kmZKulzS74rt63nervYk6X3F/97dkr4vaZK3qUTStyStk3R33bw+tyElXy6u23yXpONGoiYH39jQBXwgIo4GTgLeVVzv9Dzg+og4Cri+mLbkPcCKuunPAl8o2mojcE5Dqhp9vgT8d0Q8E3guqc28XdWRNA/4W6AtIo4hnZ3+WrxN1VwMnNFrXn/b0IuBo4phMfCvI1GQg28MiIiHI+KOYvwJ0pvTPNI1Ty8pFrsEeHljKhxdJB0K/CXwzWJawCnA5cUibitA0nTgL4CLACLiyYjYhLervjQDkyXzb4npAAAFXUlEQVQ1A1OAh/E2BUBE3AQ81mt2f9vQmcB3IrkVmClpznDX5OAbYyQtAI4FbgNaI+JhSOEIHNy4ykaVLwIfBnYX0wcBmyKiq5j29WKTI4H1wLeL3cLflDQVb1c9RMRq4ALgQVLgPQ4sw9vUQPrbhiq5drODbwyR1AJcAbw3IjY3up7RSNJLgXURsax+dh+L+ns+qRdzHPCvEXEssIXMd2v2pTg+dSZwBDAXmEraZdebt6nBVfK/6OAbIySNJ4XedyPiymL22tpuguJ2XaPqG0VOBl4maSXpJ7BOIfUAZxa7qaCf68VmaBWwKiJuK6YvJwWht6ueTgMeiIj1EbETuBJ4Pt6mBtLfNlTq2s37ysE3BhTHqC4CVkTE5+vuuhp4UzH+JuAHVdc22kTE30XEoRGxgHQCwg0R8TrgRuCVxWJuKyAiHgEekvSMYtapwG/xdtXbg8BJkqYU/4u1dvI21b/+tqGrgTcWZ3eeBDxe2yU6nHzlljFA0p8DvwD+h+7jVh8lHee7DDiM9M/5qojofZA5W5LagQ9GxEslHUnqAR4I3Am8PiJ2NLK+0UDSQtJJQBOA+4E3kz4we7uqI+ljwGtIZ1jfCbyVdGwq+21K0veBdtJPD60F/hH4L/rYhooPDl8hnQW6FXhzRCwd9pocfGZmlhPv6jQzs6w4+MzMLCsOPjMzy4qDz8zMsuLgMzOzrDj4zAYgaZek5XXDgkbXNJyKXxK4S9L7es2/WNIDks4tps+X9MFey6yUNGuA575RUqektpGp3mxomgdfxCxr2yJiYX93Smquux7jfkXSIcDzI+Lwfhb5UERc3s99g4qIRZI6hvp4s5HiHp/ZXpJ0tqT/lPRD4KfFvA9J+nXRe/pY3bJ/L+n3kn5W9K4+WMzvqPWEJM0qLqFW+53Az9U919uL+e3FY2q/jffd4su+SDpB0q8k/UbS7ZKmSfpF8eXzWh2/lPRnvV7KT4GDi57sC/ahPc6t6xE/IOnGoT6XWRXc4zMb2GRJy4vxByLiFcX484A/K642cTrp98NOJF1k92pJf0G6qPNrSb+W0QzcQbpq/0DOIV2m6QRJE4FfSvppcd+xwLNJ1y78JXCypNuBS4HXRMSvi58S2ka62srZwHslPR2YGBF39VrXy4BrBurR9vI+Sa+vm54LEBH/Bvxbcb3YG4DP9/Vgs9HCwWc2sP52dV5Xd5mu04vhzmK6hRSE04CrImIrgKSrS6zvdODPJNWu8TijeK4ngdsjYlXxXMuBBaSfwHk4In4NUPtVDkn/CfyDpA8BbyH9GOi++kJEXFCbqPVS63yJdO3THw7DusxGjIPPbGi21I0L+HREfL1+AUnvpf+fVOmi+1DDpF7P9e6I+Emv52oH6q/zuIv0/6u+1hERWyVdR/q5nFcDI3qCiaSzgcOBvxnJ9ZgNBx/jM9t3PwHeUvweIpLmSToYuAl4haTJkqYBf1X3mJXA8cX4K3s91zuK3YZIenrx46/9+R0wV9IJxfLT6n4K55vAl4Ffj+RFpCUdD3yQdBHm3YMtb9Zo7vGZ7aOI+Kmko4FbivNNOkkhcIekS4HlwJ9Iv6BRcwFwmaQ3kI6L1XyTtAvzjuLklfXAywdY95OSXgP8i6TJpON7pwGdEbFM0mbg28P0UvvzN6RfILixeP1LI+KtI7xOsyHzrzOYVUTS+aRAumCwZYdpfXOBDuCZe9sTk3Qx6cSXIX+doXieDtJPPw37T8uYDZV3dZqNQZLeSPo9xr8f4u7Hx4FP1L7APsQabgSOBHYO9TnMRoJ7fGZmlhX3+MzMLCsOPjMzy4qDz8zMsuLgMzOzrDj4zMwsK/8f/FAXsv4imqMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Optimize dimensionless, anelastic coefficients Yl\n",
    "# -------------------------------------------------\n",
    "\n",
    "# Define bound constraints for DE algorithm\n",
    "bounds = [(0.0, 1), (0.0, 1), (0.0, 1), (0.0, 1)]\n",
    "\n",
    "# Optimize Q-model by Differential Evolution\n",
    "DE_result = differential_evolution(obj_Yl, bounds)\n",
    "print('  Final obj_Qinv = ', DE_result.fun)\n",
    "\n",
    "# Calculate optimum Q model\n",
    "f, Qinv_GMB = calc_Q(DE_result.x)\n",
    "\n",
    "# Store and display optimized Yl, wl values\n",
    "Yl = np.zeros(L)\n",
    "Yl = DE_result.x\n",
    "print('Yl = ', Yl)\n",
    "print('wl = ', wl)\n",
    "\n",
    "# Calculate Q(omega)\n",
    "Q = 1 / Qinv_GMB\n",
    "\n",
    "# Define figure size\n",
    "rcParams['figure.figsize'] = 7, 5\n",
    "\n",
    "# plot stress-strain relation\n",
    "plt.plot(f, Q, 'r-',lw=3,label=\"Generalized Maxwell model\") \n",
    "plt.title(r'$Q(\\omega)$ for Generalized Maxwell model')\n",
    "plt.xlabel('Frequency f [Hz]')\n",
    "plt.ylabel(r'$Q$ []')\n",
    "plt.ylim(0,2*Qopt)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we define the modelling parameters ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "xmax = 500.0 # maximum spatial extension of the 1D model in x-direction (m)\n",
    "tmax = 0.502   # maximum recording time of the seismogram (s)\n",
    "\n",
    "vs0  = 580.   # S-wave speed in medium (m/s)\n",
    "rho0 = 1000.  # Density in medium (kg/m^3)\n",
    "\n",
    "# acquisition geometry\n",
    "xr = 330.0 # x-receiver position (m)\n",
    "\n",
    "xsrc = 250.0 # x-source position (m)\n",
    "\n",
    "f0   = 40. # dominant frequency of the source (Hz)\n",
    "t0   = 4. / f0 # source time shift (s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we first test the elastic part of the 1D SH code and compare it with the analytical solution and finally run the viscoelastic SH code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of 2D finite difference with analytical solution for homogeneous Vs model\n",
    "\n",
    "In a previous [exercise](https://danielkoehnsite.files.wordpress.com/2018/04/ex_02_tew21.pdf) you proved that the analytical solutions for the homogeneous 1D acoustic and 1D elastic SH problem, beside a density factor $1/\\rho_0$, are actual identical . In the function below we solve the homogeneous 1D SH problem by centered 2nd order spatial/temporal difference operators and compare the numerical results with the analytical solution: \n",
    "\n",
    "\\begin{equation}\n",
    "u_{y,analy}(x,t) = G_{1D} * S \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "with the 1D Green's function:\n",
    "\n",
    "\\begin{equation}\n",
    "G_{1D}(x,t) = \\dfrac{1}{2 \\rho_0 V_{s0}}H\\biggl((t-t_s)-\\dfrac{|r|}{V_{s0}}\\biggr), \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "where $H$ denotes the Heaviside function, $r = \\sqrt{(x-x_s)^2}$ the source-receiver distance (offset) and $S$ the source wavelet. Keep in mind that the stress-velocity code computes the **particle velocities** $\\mathbf{v_{y,analy}}$, while the analytical solution is expressed in terms of the **displacement** $\\mathbf{u_{y,analy}}$. Therefore, we have to take the first derivative of the analytical solution, before comparing the numerical with the analytical solution:\n",
    "\n",
    "\\begin{equation}\n",
    "v_{y,analy}(x,t) = \\frac{\\partial u_{y,analy}}{\\partial t} \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "To implement the 2D SH code, we first introduce functions to update the particle velocity $v_y$ ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Particle velocity vy update\n",
    "# ---------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_vel(vy, syx, dx, dt, nx, rho):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "            \n",
    "        # Calculate spatial derivatives            \n",
    "        syx_x = (syx[i] - syx[i - 1]) / dx\n",
    "\n",
    "        # Update particle velocities\n",
    "        vy[i] = vy[i] + (dt/rho[i]) * syx_x                \n",
    "                    \n",
    "    return vy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... update the shear stress component $\\sigma_{yx}$ for the elastic medium ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shear stress syx updates (elastic)\n",
    "# ----------------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_stress(vy, syx, dx, dt, nx, mux):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "            \n",
    "        # Calculate spatial derivatives            \n",
    "        vy_x = (vy[i + 1] - vy[i]) / dx\n",
    "\n",
    "        # Update shear stress\n",
    "        syx[i] = syx[i] + dt * mux[i] * vy_x\n",
    "                    \n",
    "    return syx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and a function for the shear stress update $\\sigma_{yx}$ for the viscoelastic case ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shear stress syx updates (viscoelastic)\n",
    "# ---------------------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_stress_visc(vy, syx, xi, Yl, wl, L, dx, dt, nx, mux):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "            \n",
    "        # Calculate spatial derivatives            \n",
    "        vy_x = (vy[i + 1] - vy[i]) / dx        \n",
    "        \n",
    "        # Calculate sum over memory variables\n",
    "        xi_sum = 0.0\n",
    "        for l in range(0, L):\n",
    "            xi_sum += Yl[l] * xi[i,l]\n",
    "        \n",
    "        # Update shear stress \n",
    "        # Note that the factor 0.5 in front of the memory variables sum \n",
    "        # is due to the arithmetic averaging of the xi variables at \n",
    "        # the staggered time steps\n",
    "        syx[i] = syx[i] + dt * mux[i] * (vy_x - 0.5 * xi_sum)\n",
    "        \n",
    "        # Update memory variables\n",
    "        xi_sum = 0.0\n",
    "        for l in range(0, L):\n",
    "            xi[i,l] += dt * wl[l] * (vy_x - xi[i,l])\n",
    "            \n",
    "            # After xi update calculate new xi_sum ...\n",
    "            xi_sum += Yl[l] * xi[i,l]\n",
    "        \n",
    "        # ... and finish stress update\n",
    "        syx[i] = syx[i] - dt * mux[i] * xi_sum\n",
    "        \n",
    "    return syx, xi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and harmonically averaging the (unrelaxed) shear modulus ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Harmonic averages of shear modulus\n",
    "# ----------------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def shear_avg(mu, nx, mux):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "            \n",
    "        # Calculate harmonic averages of shear moduli        \n",
    "        mux[i] = 2 / (1 / mu[i + 1] + 1 / mu[i])\n",
    "            \n",
    "    return mux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can assemble the main FD code ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "code_folding": [
     42
    ]
   },
   "outputs": [],
   "source": [
    "# 2D SH viscoelastic wave propagation (Finite Difference Solution) \n",
    "# ----------------------------------------------------------------\n",
    "def FD_1D_visc_SH_JIT(dt,dx,f0,xsrc,Yl,wl,L,mode):\n",
    "        \n",
    "    nx = (int)(xmax/dx) # number of grid points in x-direction\n",
    "    print('nx = ',nx)    \n",
    "            \n",
    "    nt = (int)(tmax/dt) # maximum number of time steps            \n",
    "    print('nt = ',nt)\n",
    "    \n",
    "    ir = (int)(xr/dx)      # receiver location in grid in x-direction    \n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "\n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of a Gaussian\n",
    "    src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "    # Analytical solution\n",
    "    # -------------------\n",
    "    G        = time * 0.\n",
    "    vy_analy = time * 0.\n",
    "\n",
    "    # Initialize coordinates\n",
    "    # ----------------------\n",
    "    x    = np.arange(nx)\n",
    "    x    = x * dx       # coordinates in x-direction (m)\n",
    "    \n",
    "    # calculate source-receiver distance\n",
    "    r = np.sqrt((x[ir] - x[isrc])**2)\n",
    "    \n",
    "    for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "        if (time[it] - r / vs0) >= 0:\n",
    "            G[it] = 1. / (2 * rho0 * vs0)\n",
    "    Gc   = np.convolve(G, src * dt)\n",
    "    Gc   = Gc[0:nt]\n",
    "    \n",
    "    # compute vy_analy from uy_analy\n",
    "    for i in range(1, nt - 1):\n",
    "        vy_analy[i] = (Gc[i+1] - Gc[i-1]) / (2.0 * dt)           \n",
    "    \n",
    "    # Initialize empty wavefield arrays\n",
    "    # ---------------------------------\n",
    "    vy    = np.zeros(nx) # particle velocity vy\n",
    "    syx   = np.zeros(nx) # shear stress syx     \n",
    "\n",
    "    # Initialize model (assume homogeneous model)\n",
    "    # -------------------------------------------    \n",
    "    vs    = np.zeros(nx)\n",
    "    vs    = vs + vs0      # initialize wave velocity in model\n",
    "    \n",
    "    rho   = np.zeros(nx)\n",
    "    rho   = rho + rho0    # initialize wave velocity in model\n",
    "    \n",
    "    # calculate shear modulus\n",
    "    # -----------------------\n",
    "    mu    = np.zeros(nx)\n",
    "    mu    = rho * vs ** 2 \n",
    "    \n",
    "    # Estimate unrelaxed shear modulus in viscoelastic case\n",
    "    # -----------------------------------------------------\n",
    "    if(mode=='visc'):\n",
    "        mu1 = mu / (1 - np.sum(Yl))\n",
    "    \n",
    "    # harmonic average of shear moduli\n",
    "    # --------------------------------\n",
    "    if(mode=='elast'):\n",
    "        mux   = mu # initialize harmonic average mux\n",
    "    \n",
    "    if(mode=='visc'):\n",
    "        mux   = mu1 # initialize harmonic average mux \n",
    "    \n",
    "    mux = shear_avg(mu, nx, mux)\n",
    "    \n",
    "    # Initialize memory variables\n",
    "    # ---------------------------\n",
    "    if(mode=='visc'):\n",
    "        xi = np.zeros((nx,L))\n",
    "    \n",
    "    # Initialize empty seismogram\n",
    "    # ---------------------------\n",
    "    seis = np.zeros(nt) \n",
    "    \n",
    "    # Time looping\n",
    "    # ------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # Update particle velocity vy\n",
    "        # ---------------------------\n",
    "        vy = update_vel(vy, syx, dx, dt, nx, rho)\n",
    "\n",
    "        # Add Source Term at isrc\n",
    "        # ------------------------------\n",
    "        # Absolute particle velocity w.r.t analytical solution\n",
    "        vy[isrc] = vy[isrc] + (dt * src[it] / (rho[isrc] * dx))\n",
    "        \n",
    "        # Update shear stress syx, syz\n",
    "        # ----------------------------\n",
    "        if(mode=='elast'):\n",
    "            syx = update_stress(vy, syx, dx, dt, nx, mux)                \n",
    "            \n",
    "        if(mode=='visc'):\n",
    "            syx, xi = update_stress_visc(vy, syx, xi, Yl, wl, L, dx, dt, nx, mux)                           \n",
    "            \n",
    "        # Output of Seismogram\n",
    "        # -----------------\n",
    "        seis[it] = vy[ir]   \n",
    "        \n",
    "    # Compare FD Seismogram with analytical solution\n",
    "    # ---------------------------------------------- \n",
    "    # Define figure size\n",
    "    rcParams['figure.figsize'] = 12, 5\n",
    "    \n",
    "    if(mode=='elast'):\n",
    "        label = \"Elastic FD solution\"\n",
    "        \n",
    "    if(mode=='visc'):\n",
    "        label = \"Viscoelastic FD solution (Q = \" + str(Qopt) + \")\"    \n",
    "    \n",
    "    plt.plot(time, seis, 'b-',lw=3,label=label) # plot FD seismogram\n",
    "    Analy_seis = plt.plot(time,vy_analy,'r--',lw=3,label=\"Elastic analytical solution\") # plot analytical solution\n",
    "    plt.xlim(time[0], time[-1])\n",
    "    plt.title('Seismogram')\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... run the elastic FD code and compare the result with the analytical solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nx =  500\n",
      "nt =  502\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvcAAAFNCAYAAACExkA1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl4VdXZ9/HvnTmEkJAwEyAoyDwPDqigqKBQEUHBqepTpVrFPo+tFdu+TtW2WhyrbZ3Rlio4IVWpI3FCZRAEmWSGACKEAAmQkGG9f5zNyRwOZNhJ+H2uKxdrr73XWvcJm8N91ll7b3POISIiIiIi9V+Y3wGIiIiIiEj1UHIvIiIiItJAKLkXEREREWkglNyLiIiIiDQQSu5FRERERBoIJfciIiIiIg2EknsRESnBzLLN7AS/4xARkaOn5F5EpIEys9PNbJ6Z7TWz3Wb2hZkNOlI751xj59z62ohRRESqV4TfAYiISPUzsybA28CNwEwgCjgDyPUzrupgZuHOuQK/4xARqYs0cy8i0jCdBOCce9k5V+CcO+ice985txTAzP7HzFaaWaaZvWdmHQ43NDNnZp288gVmtsLMssxsq5n92qsfZmbpZvYbM/vRzLab2UXe8d973xT8tlif0Wb2qJlt834eNbPoYvt/4/WxzcyuKxXDNDP7u5m9a2b7gbPMbJSZLTazfWa2xczuLtZXqtf+Wm9fppndYGaDzGypme0xsydq9tcvIuIPJfciIg3T90CBmb1oZuebWdPDO8zsIuC3wMVAc+Az4OUK+nkO+LlzLh7oCXxcbF8rIAZoC9wJPANcCQwg8C3BncXW7v8OOAXoC/QBBgO/9+IZCdwKnAN0AoaWE8flwP1APPA5sB/4KZAIjAJu9F5XcScDnYEJwKNeDOcAPYBLzay8cURE6jUl9yIiDZBzbh9wOuAIJN07zWy2mbUEfg78yTm30jmXD/wR6Ft89r6YPKC7mTVxzmU6574pte9+51we8ArQDHjMOZflnFsOLAd6e8deAdzrnPvRObcTuAe4ytt3KfCCc265c+6At6+0t5xzXzjnCp1zOc65NOfcMm97KYEPJ6WT9T94x75P4MPAy974Wwl8oOkX2m9TRKT+UHIvItJAecn7Nc65FAKz7m0IzGB3AB7zlqfsAXYDRmAGvrRxwAXAJjP7xMxOLbYvo9ja94PenzuK7T8INPbKbYBNxfZt8uoO79tSbF/xcrl1Znaymc01s51mthe4gcCHi+JKx1JRbCIiDYaSexGR44BzbhUwjUCSv4XAUpvEYj+xzrl55bRb4JwbA7QAZhG4OPdYbCPwoeKw9l4dwHYgpdi+duW9hFLb/wZmA+2ccwnAPwh8QBEROa4puRcRaYDMrKuZ/crMUrztdsBlwFcEEuE7zKyHty/BzC4pp48oM7vCzBK8pTf7gGO9S83LwO/NrLmZNSOwRv9f3r6ZwLVm1s3MGnn7jiQe2O2cyzGzwQTW5IuIHPeU3IuINExZBC4o/dq7w8xXwHfAr5xzbwIPAK+Y2T6v/vwK+rkK2OgddwOBC2aPxX3AQmApsAz4xqvDOTcHeByYC6wFvvTaVHbbzl8A95pZFoEPA8f6jYKISINizpX+plNERMQ/ZtaNwAeOaO+CXxERCZFm7kVExHdmNtZbBtSUwLcK/1FiLyJy9JTci4hIXfBzYCewjsC6/hv9DUdEpH7SshwRERERkQZCM/ciIiIiIg2EknsRERERkQYiwu8A6qPExETXqVMnv8OQOmb//v3ExcX5HYbUITonpDw6L6Q0nRNSWvFzYtGiRbucc81Dbavk/hi0bNmShQsX+h2G1DFpaWkMGzbM7zCkDtE5IeXReSGl6ZyQ0oqfE2a26WjaalmOiIiIiEgDoeReRERERKSBUHIvIiIiItJAaM29iIiIHNfy8vJIT08nJyfHl/ETEhJYuXKlL2NL3RETE0NKSgqRkZFV6kfJvYiIiBzX0tPTiY+PJzU1FTOr9fGzsrKIj4+v9XGl7nDOkZGRQXp6Oh07dqxSX74uyzGzkWa22szWmtmUcvZHm9kMb//XZpZabN8dXv1qMxtxpD7NrKPXxxqvzyiv/gYzW2ZmS8zsczPrXrOvWkREROqSnJwckpOTfUnsRQDMjOTk5Gr59si35N7MwoEngfOB7sBl5STWPwMynXOdgEeAB7y23YGJQA9gJPA3Mws/Qp8PAI845zoDmV7fAP92zvVyzvUFHgQerpEXLCIiInWWEnvxW3Wdg37O3A8G1jrn1jvnDgGvAGNKHTMGeNErvwYMt8ArHwO84pzLdc5tANZ6/ZXbp9fmbK8PvD4vAnDO7Ss2Xhzgqvl1ioiIiFQqPDycvn37Bn/+/Oc/AzBs2LBjerbOrFmzWLFiRXD7zjvv5MMPPwyp7caNG4mNjS0Rz6FDh5g2bRrNmzenX79+dO7cmREjRjBv3ryjjq20tLQ0Ro8eXekxe/bs4W9/+1twe9u2bYwfP77KYzdEfq65bwtsKbadDpxc0THOuXwz2wske/VflWrb1iuX12cysMc5l1/O8ZjZTcCtQBSBDwFlmNkkYBJA8+bNSUtLC+U1ynEkOztb54WUUCvnhHNEZGVBbh75zZJAs491nt4r6p6EhASysrJ8G7+goIDY2Fg+++yzEvVZWVkUFBSwf//+o47v1VdfZeTIkbRr1w6A2267LdjnkWRnZ9OxY8cS8eTm5pKTk8PYsWN56KGHAPj0008ZO3Ys77zzDl26dDmq+Io7cOAA+fn5lcaWnp7OE088wVVXXQVAfHw8L7zwgq9/bzUhJyeHtLS0Kr1P+Jncl/c/UOlZ84qOqai+vG8iKjs+UHDuSeBJM7sc+D1wdZmDnXsaeBqgS5cuTk+Sk9L0hEEprabPiR0z0sidNJn2+5YD8F3i6ex7+1NOG6IEvy7Te0Xds3LlSl8vaD2coJYXQ3h4OHFxccTHx3PjjTeyYMECDh48yPjx47nnnnsAmDJlCrNnzyYiIoLzzjuPiy++mDlz5jBv3jweeughXn/9df7whz8wevRoxo8fz4IFC/jlL3/J/v37iY6O5qOPPioxduPGjQkLCysTT0xMDFFRUcH6UaNG8fOf/5zp06fzyCOPlDj21Vdf5Z577iE8PJyEhAQ+/fRTcnJyuPHGG1m4cCERERE8/PDDnHXWWTRq1IiIiAji4+O5++67ady4Mb/+9a8B6NmzJ2+//Tb33XcfGzZs4IwzzuDcc8/lpptuYvTo0Xz33XcV9jtt2jRmz57NgQMHWLduHWPHjuXBBx+svr+4GhATE0O/fv2q9D7hZ3KfDrQrtp0CbKvgmHQziwASgN1HaFte/S4g0cwivNn78saCwDKevx/TqxERqUWrb3mSLn+9uUTd7D1n8P/ONGbMAH1bLVK/HDx4kL59+wa377jjDiZMmFDimPvvv5+kpCQKCgoYPnw4S5cuJSUlhTfffJNVq1ZhZuzZs4fExEQuvPDCYDJf3KFDh5gwYQIzZsxg0KBB7Nu3j9jY2DLxrFu3LhjPkCFDePLJJ8uNu3///jz11FNl6u+9917ee+892rZty549ewCCfSxbtoxVq1Zx3nnn8f3334f0+/nzn//Md999x5IlS4DA0qHDKut3yZIlLF68mOjoaLp06cLkyZOD32Y0VH4m9wuAzmbWEdhK4ALZy0sdM5vALPqXwHjgY+ecM7PZwL/N7GGgDdAZmE9ghr5Mn16buV4fr3h9vgVgZp2dc2u88UYBaxARqcPSn3qHTn+9pURdFo2Zyq8pLIQrroBWreD0030KUKQeq8mVba6Sq/piY2ODiWtFZs6cydNPP01+fj7bt29nxYoVdO/enZiYGK677jpGjRp1xLXrq1evpnXr1gwaNAiAJk2alHvciSeeeMR4IHALx/IMGTKEa665hksvvZSLL74YgM8//5zJkycD0LVrVzp06BBycl+ZyvodPnw4CQkJAHTv3p1NmzY1+OTetwtqvRn0m4H3gJXATOfccjO718wu9A57Dkg2s7UE1sRP8douB2YCK4D/Ajc55woq6tPr63bgVq+vZK9vgJvNbLmZLfHGKLMkR0Skrsjd8iPxv7iScAoBWBw5mC9mbmXbt7toflISAIcOwfXXQ16en5GKSHXasGEDU6dO5aOPPmLp0qWMGjWKnJwcIiIimD9/PuPGjWPWrFmMHDmy0n6cc9V6Z6DFixfTrVu3MvX/+Mc/uO+++9iyZQt9+/YlIyOjwg8CxUVERFBYWBjcDuXWkJX1Gx0dHSyHh4eTn59f4bENha/3uXfOveucO8k5d6Jz7n6v7k7n3GyvnOOcu8Q518k5N9g5t75Y2/u9dl2cc3Mq69OrX+/10cnrM9er/6Vzrodzrq9z7qxiHwZEROqc1WOnkFAY+Ip7s7UnYs5/GHJJG7r0jua99+DwEtlVq+DZx/b7GKmIVKd9+/YRFxdHQkICO3bsYM6cQOqTnZ3N3r17ueCCC3j00UeDs+3x8fHlXmzatWtXtm3bxoIFC4DAev9jTXg/+eQTnn76aa6//voy+9atW8fJJ5/MvffeS7NmzdiyZQtnnnkm06dPB+D7779n8+bNZS7ETU1N5ZtvvgHgm2++YcOGDZW+HiCkfo8nekKtiEg9sXPOQnoveiG4/c31/+Ci4S2C26mp8Pvfw0O37+Be7mTU7XM4OGkdsU2q9ihzkeNJCJPLNaL0mvuRI0cGb4cJ0KdPH/r160ePHj044YQTGDJkCBBIzseMGUNOTg7OueCFrRMnTuT666/n8ccf57XXXgv2ExUVxYwZM5g8eTIHDx4kNjaWDz/8kMaNG4cU54wZM/j88885cOAAHTt25PXXXy935v62225jzZo1OOcYPnw4ffr0oWvXrtxwww306tWLiIgIpk2bVmJmHWDcuHG89NJL9O3bl0GDBnHSSScBkJyczJAhQ+jZsyfnn38+N910U7DNL37xiyP2ezyxUL4ikZK6dOniVq9e7XcYUsfoDhhSWnWfE8s6jaXXulkAfBQ/hjMzZhFZKm/PzXFsje/CCfmBy4fSbnqVYU/o6tq6RO8Vdc/KlSvLTVBrS1ZWlq9365G64/C5WPx9wswWOecGhtqHr8tyREQkNJnfbaX7utnB7fA/3lcmsQeIjjG2D70suB37z7J3sRARkYZLyb2ISD3w1zfa0pPvmMbVvN9kHENv6lnhsV2nXkeB9/Z+8r4PWf7W2toKU0REfKbkXkSkjsvPh6eeglV041qmsfOJmZXeri+5bzuWtB0V3N7y0MxaiFJEROoCJfciInXcu+/CNu+xey1bwqUTj/zWbRMnBstt57/h20WCIiJSu5Tci4jUcU8/XVS+9lrKXWtfWvfbRnGIwIG9chex8r+baig6ERGpS5Tci4jUYTs/XUm7d58iiQwArrsutHYxLRNY0Xp4cHvTo2/WRHgiIlLHKLkXEanDNt/9PH93N/ADrXjyhIc48cTQ2xZceHGwnPTlOzUQnYhUl/DwcPr27Rv8OXyP+2HDhrFw4cKj7m/WrFmsWLEiuH3nnXfy4YcfVlu8obr77ruZOnVqpcc8+uijHDhwILh9wQUXsGfPnqMea9q0adx8881H3a64jRs30rNnxTcsOOyPf/xjie3TTjutSuNWJyX3IiJ1VWEh7T5/GYBI8uk4qvtRNe90c9Fj6Htnfc7ubUd+jLuI+CM2NpYlS5YEf6ZMmVKl/kon9/feey/nnHNOVcOsEaWT+3fffZfExEQfIzqy0sn9vHnzfIqkLCX3IiJ1VPrLn9EibysAO2nGyb87uv+YE3q2Y3NMZwC20I4Fb2yp9hhFpPbceOONDBw4kB49enDXXXcF66dMmUL37t3p3bs3v/71r5k3bx6zZ8/mtttuo2/fvqxbt45rrrkm+KTaBQsWcNppp9GnTx8GDx5MVlZWiXGys7MZPnw4/fv3p1evXrz11ltAYFa7W7duXH/99fTo0YPzzjuPgwcPAvDMM88waNAg+vTpw7hx40ok6wDr1q2jf//+we01a9YwYMAAHn/8cbZt28ZZZ53FWWedBUBqaiq7du0C4KWXXqJ379706dOHq666CoD//Oc/nHzyyfTr149zzjmHHTt2VPp7++STT4LfiPTr14+srCycc9x222307NmTXr16MWPGjDLtSn8TMHr0aNLS0pgyZUrwicJXXHEFQPAJvxX1e/ihVOPHj6dr165cccUV1NSDZCNqpFcREamy7X97kxSvvCD1Ei5oGcKVtKXMmfAi97/Yli2054blMKJ6QxSRanI4WTzsjjvuYMKECSWOuf/++0lKSqKgoIDhw4ezdOlSUlJSePPNN1m1ahVmxp49e0hMTOTCCy9k9OjRjB9f8gnVhw4dYsKECcyYMYNBgwaxb98+YmNjSxwTExPDm2++SZMmTdi1axennHIKF154IRBIyl9++WWeeeYZLr30Ul5//XWuvPJKLr74Yq6//noAfv/73/Pcc88xefLkYJ8nnngiCQkJLFmyhL59+/LCCy9wzTXXMHnyZB5++GHmzp1Ls2bNSsSxfPly7r//fr744guaNWvG7t27ATj99NP56quvMDOeffZZHnzwQR566KEKf7dTp07lySefZMiQIWRnZxMTE8Mbb7zBkiVL+Pbbb9m1axeDBg3izDPPDOnv6s9//jNPPPEES5YsKbOvsn4XL17M8uXLadOmDUOGDOGLL77g9NNPD2nMo6GZexGRusg52i56K7gZdenYY+qmyzWnsoX2AHzwQbVEJtLw3X03mIX2M2lS2faTJpU85u67jzhk6WU5pRN7gJkzZ9K/f3/69evH8uXLWbFiBU2aNCEmJobrrruON954g0aNGlU6zurVq2ndujWDBg0CoEmTJkRElJzrdc7x29/+lt69e3POOeewdevW4Ox4x44dgx9CBgwYwMaNGwH47rvvOOOMM+jVqxfTp09n+fLlZca+7rrreOGFFygoKGDGjBlcfvnllcb68ccfM378+GDSn5SUBEB6ejojRoygV69e/OUvfyl3rOKGDBnCrbfeyuOPP86ePXuIiIjg888/57LLLiM8PJyWLVsydOhQFixYUGk/oais38GDB5OSkkJYWBh9+/YN/u6qm5J7EZE6aNfcZbTJ3QjAXprQ//+GHlM/p54Khyfl1q0rul++iNQvGzZsYOrUqXz00UcsXbqUUaNGkZOTQ0REBPPnz2fcuHHMmjWLkSNHVtqPcw6r7Cl4wPTp09m5cyeLFi1iyZIltGzZkpycwDU70dHRwePCw8PJz88H4JprruGJJ55g2bJl3HXXXcHjixs3bhxz5szh7bffZsCAASQnJx9TrJMnT+bmm29m2bJlPPXUU+WOVdyUKVN49tlnOXjwIKeccgqrVq0KaUlMREQEhYWFwe0jjXM45opU9LurbkruRUTqoI2PFc3aL2xxAUmtoo6pn+ho8CboAPjyy6pGJiJ+2LdvH3FxcSQkJLBjxw7mzJkDBNbH7927lwsuuIBHH300uFQkPj6+zFp6gK5du7Jt27bgbHJWVlaZJHPv3r20aNGCyMhI5s6dy6ZNR35ORlZWFq1btyYvL4/p06eXe0xMTAwjRozgxhtv5Nprrw3WVxTr8OHDmTlzJhkZgVsBH16Ws3fvXtq2bQvAiy++eMTY1q1bR69evbj99tsZOHAgq1at4swzz2TGjBkUFBSwc+dOPv30UwYPHlyiXWpqKkuWLKGwsJAtW7Ywf/784L7IyEjy8vLKjBVKvzVNa+5FROqghE+KkvvcEWOq1NcZJx/i4KdLOJUv2fdiexh3bEt8RI4bd98d0lKaCj39dMmnz4Wg9Jr7kSNHBm+HCdCnTx/69etHjx49OOGEExgyZAgQSKrHjBlDTk4OzjkeeeQRACZOnMj111/P448/HryQFiAqKooZM2YwefJkDh48SGxsLB9++GHwglCAK664gp/85CcMHDiQvn370rVr1yPG/4c//IGTTz6ZDh060KtXr3KT9cN9v/HGG5x33nnBukmTJnH++efTunVr5s6dG6zv0aMHv/vd7xg6dCjh4eH069ePadOmcffdd3PJJZfQtm1bTjnlFDZs2FBpbI8++ihz584lPDyc7t27c/755xMVFcWXX35Jnz59MDMefPBBWrVqVWKpzJAhQ+jYsSO9evWiZ8+eJS4InjRpEr1796Z///4lPsyMHTu23H5XrVp1xN9hdbGaulK3IevSpYtbvXq132FIHXP4SniRw471nNi/Op24ru0AOEQk25bsJLVPwjHHsfh/X6TfY9cA8EXiKIZkvn3MfUnV6b2i7lm5ciXdunXzbfysrCzi4+N9G782TZ06lb179/KHP/zB71DqpMPnYvH3CTNb5JwbGGofmrkXEaljvp86m35eeWHjYZxWhcQeoMNlp8JjgXLXPV+Sc9ARE1v5mlsRkeo2duxY1q1bx8cff+x3KA2aknsRkTrmyeyryaA1Y3iL+LPPrnJ/SYM7kxmeTNOCDJLZzaI3v2fA5V2qIVIRkdC9+eabfodwXNAFtSIidUhBAbz5fhyzGMu1TKPD//tp1Ts1Y1PrU4Obu975uup9iohInaTkXkSkDlm4ELwbQtCmDQwYUD395vYuumWO++ab6ulUpAHRNYjit+o6B5Xci4jUIcUfNHXuuYHn31SHJsOK7vLQbJOSe5HiYmJiyMjIUIIvvnHOkZGRQUxMTJX70pp7EZE6ZMsbC4igL/lEcu651ddvh7H94TeBcpeDi9mfVUhcvOZ3RABSUlJIT09n586dvoyfk5NTLUmd1G8xMTGkpKRUuR8l9yIidcT+tdt5avFgptKY9zmPIWe/BlTP1H2jE1uzK7wlzQp2EE8237y9hv6X6aJaEQg8kKhjx46+jZ+Wlka/fv2OfKBICDRtIyJSR6z5x4cAxJNNu7hMWrWuxttVmrG1ZdHSnB//q6U5IiINkZJ7EZE6Iu/dD4PlXX3Pqfb+c3sUJfdukZJ7EZGGSMtyRETqAudov6YouU8cX/3JfcyIobz5wXKW0YvNh87h/GofQURE/KbkXkSkDtj52Spa5m8DIJNE+vxPNd0Ds5iUa8+lz68DV+lGb4Z/5EOE/hcQEWlQtCxHRKQO2PhM0T0wlzU7m7gm4dU+RlIStG0bKOfmwtq11T6EiIj4TMm9iEgdEJFWtCRn/6nVvyTnsF69ispLl9bYMCIi4hMl9yIiPnOH8ui0NS243eantZPcL1tWY8OIiIhPlNyLiPhs46sLiHdZAGyx9vS8qFONjTU0/hue51oWMJCTp99SY+OIiIg/fE3uzWykma02s7VmNqWc/dFmNsPb/7WZpRbbd4dXv9rMRhypTzPr6PWxxuszyqu/1cxWmNlSM/vIzDrU7KsWESlpx7+K1tuvbn8u4RHVeH/7Uro03821TGMgi2i7bX6NjSMiIv7wLbk3s3DgSeB8oDtwmZl1L3XYz4BM51wn4BHgAa9td2Ai0AMYCfzNzMKP0OcDwCPOuc5Aptc3wGJgoHOuN/Aa8GBNvF4RkYp8sqMr73EeB4nBDa+5JTkA7S4oWpfTOfc79mcV1uh4IiJSu/ycuR8MrHXOrXfOHQJeAcaUOmYM8KJXfg0Ybmbm1b/inMt1zm0A1nr9ldun1+Zsrw+8Pi8CcM7Ndc4d8Oq/AlJq4LWKiJQrLw/+8P0ERvIeTcnkxF9dVKPjRbdvSWZ4MgCN2c/6T9NrdDwREaldfib3bYEtxbbTvbpyj3HO5QN7geRK2lZUnwzs8fqoaCwIzObPOYbXIiJyTL76CvbvD5Rbp8bQsVtMjY/5Q0KXYPnHz1bX+HgiIlJ7/Hx8SXmLSl2Ix1RUX96HlcqOLxrI7EpgIDC0nGMxs0nAJIDmzZuTlpZW3mFyHMvOztZ5ISWEck48/3wqkApAjx7b+OST72s6LEhsS7fdgWL6x1+SlhZZ82NKkN4rpDSdE1JaVc4JP5P7dKBdse0UYFsFx6SbWQSQAOw+Qtvy6ncBiWYW4c3elxjLzM4BfgcMdc7llhesc+5p4GmALl26uGHDhoX8QuX4kJaWhs4LKS6Uc+K3vy0qX311G4YNa1OzQQGLBn4N618FoHnGTp23tUzvFVKazgkprSrnhJ/LchYAnb272EQRuEB2dqljZgNXe+XxwMfOOefVT/TuptMR6AzMr6hPr81crw+8Pt8CMLN+wFPAhc65H2votYqIlLFv1Tae/7Irf+VmRvEOZ59dO+M2GVi0LCdhh5bliIg0JL7N3Dvn8s3sZuA9IBx43jm33MzuBRY652YDzwH/NLO1BGbsJ3ptl5vZTGAFkA/c5JwrACivT2/I24FXzOw+AnfIec6r/wvQGHg1cN0tm51zF9bwyxcRYc0/PmIAq+nKak6OX0Fy8qhaGbfN2V2D5ZT9q8nLg0itzBERaRD8XJaDc+5d4N1SdXcWK+cAl1TQ9n7g/lD69OrXE7ibTun6mr3vnIhIBfL/+2GwvLtf7b0VxfU6gXzCiaCADmxm1dL9dB0QV2vji4hIzdETakVE/OAcqWuLkvukS2txniEqih8anRDc3Jq2pvbGFhGRGqXkXkTEB9s+WknLgsB1/Zkk0uuaAbU6/n+HPcAFvMOJrOWbvF5HbiAiIvWCknsRER9sfr5o1n55i7OJiQuv1fEPjBjLHC5gPSfy/braHVtERGqOknsRER9EflKU3B8cUvuX/nTuXFReu7bWhxcRkRqi5F5EpJYV5uZx0va04Hbbq2s/ue/Uqai8RkvuRUQaDCX3IiK1bN3L84l3WQBsDutA19GdjtCi+qWmQng4RHKIJltXcGB/6QeEi4hIfaTkXkSklv3476IlOWtTzyEs3Go9hshIWBQxmIPEsoIebPqy9APCRUSkPlJyLyJSy3K+K1oHY+f496iN6JgwwikEYNc8PalWRKQhUHIvIlKLcnJgdOa/SGEL1/ACJ910rm+x7G3VJVg+uFjJvYhIQ6DkXkSkFn3xRSDB30oKX3W5hra9k32LpaBTUXIftlbJvYhIQ6DkXkSkFn1YtNweH1fkABDTpyi5b7ztex8jERGR6qLkXkSkFn3wQVH5XP9W5ACQfHLRXXqa7VvvYyQiIlJdlNyLiNSS3St3cPaiv9CXxUSEFTJsmL/xtBnSMVhul78KBh42AAAgAElEQVSBA9mFPkYjIiLVQcm9iEgtWfPEezzIb1hMfz5KuJiEBH/jiUxuwu6wZgBEc4jNX+l2mCIi9Z2SexGRWuL++16wfKjvYB8jKfJj/AlF5a+0NEdEpL5Tci8iUgtcQSGdN74f3G5x1QgfoymS3eLEYHn/0nU+RiIiItVByb2ISC1Y//pikgt3AbDTmtP9in4+RxRQ2DEwc59BEru3HfQ5GhERqSol9yIitWDbtKIlOatSziUiqm68/e79n1tJJJNmZPBc1C/8DkdERKoowu8ARESOB02+KlqSU3hu3ViSA5DaP4m9XnnNGl9DERGRalA3po5ERBqwAzuy6J75RXC7803n+RhNSampEB4eKKenw0GtzBERqdeU3IuI1LCVf5tLJPkArIruQ5v+rXyOqEhkJLRrV7S9aZN/sYiISNUpuRcRqWEHZhWtt9/Wq+4syTmsX5sdDOFzruIltn+z3e9wRESkCpTci4jUsPYri5L7JpfUveT+7k3X8jln8BJXk/PJ136HIyIiVaALakVEatCG9Y7L8/7JebzP2WFpDJ40xO+QyshtewJsDZQLvte97kVE6jMl9yIiNeg/bxtfcSpfcSoLR97FO4l+R1SWnXgCzA+UI7foKbUiIvWZluWIiNSg2bOLyj/5iX9xVKZRzxOC5fidSu5FROozJfciIjUkOzucTz4p2h492r9YKpM8qCi5b7Ffyb2ISH2m5F5EpIasfX8/zfIDd5/p1w9SUnwOqALNB3cMltsXbCB7b4GP0YiISFUouRcRqSE9X3uW7bThawZzc69PjtzAJ2EJ8ewKbwFAFHmkf73V54hERORYKbkXEakB+bkF9P9hLgCDWcAp5zT2OaLK/di4aGnO7oVamiMiUl8puRcRqQHLn/2SJLcbgB/CWtP1sn4+R1S5rOZFyf2B75Tci4jUV0ruRURqwO6X3g6Wvz9pNGERdfvttiAltai8fpN/gYiISJXU7f9tRETqIVfo6LD4zeB2zLg6epucYsK6d2UZPXmbUazM6+R3OCIicox8Te7NbKSZrTaztWY2pZz90WY2w9v/tZmlFtt3h1e/2sxGHKlPM+vo9bHG6zPKqz/TzL4xs3wzG1+zr1hEjgerX/+OE/K+ByCbOHrdeq7PER2Z/fQqerOMn/A20wqu8jscERE5Rr4l92YWDjwJnA90By4zs+6lDvsZkOmc6wQ8Ajzgte0OTAR6ACOBv5lZ+BH6fAB4xDnXGcj0+gbYDFwD/LsmXqeIHH+2PfZqsLys/Whik2J9jCY0qalF5Q0bwDnfQhERkSrwc+Z+MLDWObfeOXcIeAUYU+qYMcCLXvk1YLiZmVf/inMu1zm3AVjr9Vdun16bs70+8Pq8CMA5t9E5txQorKkXKiLHD+eg/YLXgtthl9aPLwRbtIBGjQLlfftgzx5/4xERkWPjZ3LfFthSbDvdqyv3GOdcPrAXSK6kbUX1ycAer4+KxhIRqbLVbyyn06GVAOynEb2nXOBzRKExKzt7LyIi9U+Ej2NbOXWlvwiu6JiK6sv7sFLZ8SEzs0nAJIDmzZuTlpZ2NM3lOJCdna3zQth9zyy6euX5zc7Gls33NZ6jcXbhbkazhg5s4utnz2HfpUl+h9Qg6b1CStM5IaVV5ZzwM7lPB9oV204BtlVwTLqZRQAJwO4jtC2vfheQaGYR3ux9eWNVyjn3NPA0QJcuXdywYcOOprkcB9LS0tB5cXxzDibsiWEHuVzMG+w+eyjj6tE5EZ8/gQHMBGBO5ikMG3axzxE1THqvkNJ0TkhpVTkn/FyWswDo7N3FJorABbKzSx0zG7jaK48HPnbOOa9+onc3nY5AZ2B+RX16beZ6feD1+VYNvjYROQ4tWwavbjmFX/B3TorbRvxP+/sd0lEpaNuhqKx73YuI1Eu+JffeDPrNwHvASmCmc265md1rZhd6hz0HJJvZWuBWYIrXdjkwE1gB/Be4yTlXUFGfXl+3A7d6fSV7fWNmg8wsHbgEeMrMDh8vInJUpk8vKl/wk3Ci4urXo0QiO6cWlbcruRcRqY/8XJaDc+5d4N1SdXcWK+cQSLrLa3s/cH8ofXr16wncTad0/QICy3RERI5Zfj78859F21de6V8sx6pxj6KZ+4TdG/0LREREjln9mlYSEamjPvv3FnK27wagZUsYMeIIDeqgpP6pwXLzg5q5FxGpj5Tci4hUg7A7f892WjOTS/jNyKVE+Pq96LFJ6lc0c59SuJm9mXr8h4hIfaPkXkSkinZvymLgpteI5hCX8Bpjzj/kd0jHxOIbkxmWDEA0h9j2zQ8+RyQiIkdLyb2ISBUt/v3rxHEAgHWxPTjx0gE+R3TsdsYVzd7vXqylOSIi9Y2SexGRKkp8a1qwvGPkNYHHvdZTWUmpwfKBFRt9i0NERI6NknsRkSr47vXVDMj6BIB8wul23xU+R1Q1eW2KZu7zda97EZF6px5e8iUiUnds/91f6emVv207igHdW/saT1XtHziUp77czyY6EB59Nuf7HZCIiBwVJfciIsdox+o9nLp6WnA77re/9C+YamIXjeGGv44BYMh+n4MREZGjpmU5IiLHaMktz9OYQAa8LrYnXW88y+eIqq59+6Ly5s3+xSEiIsdGyb2IyDHIPVBAtw//GtzOvOqWen0h7WHt2hWVt24NPHlXRETqDyX3IiLH4LP/9z7tCzcCkBmWRJ8H6/eFtIdFR0OrVoFyYWEgwRcRkfpDyb2IyFEqKID/nTOCkcxhDiNZeebPiUxo5HdY1eZ3EQ8wizEspi8ZHy3xOxwRETkKuqBWROQovfIKLF8ZxnJGMi9+JBtfLfQ7pGp1asFnDOAdAD5duh7o629AIiISMs3ci4gchfx8uOeeou3//V9Iataw3kpzWxRdVXtoje51LyJSnxzxfyQza2Rm/8/MnvG2O5vZ6JoPTUSk7vnnP2HNmkA5MRFuvdXfeGqCa1/0ICvbolvmiIjUJ6FMN70A5AKnetvpwH01FpGISB11KKeQtjddxP/wHBHk8etfBxL8hiaqc9HMfcyPmrkXEalPQknuT3TOPQjkATjnDgL1/35vIiJH6aMbXuW8g2/xHNexOHwgt9xU4HdINSKhV9HMfeJezdyLiNQnoST3h8wsFnAAZnYigZl8EZHjxo8b9tPzpduC2wfOGEl8YriPEdWc5gOLkvuWOZtwzsdgRETkqISS3N8F/BdoZ2bTgY+A39RoVCIidczXF/2Jdm4LABnhzenzyh0+R1RzEru2Is+7mVozdrFn636fIxIRkVAdMbl3zn0AXAxcA7wMDHTOpdVsWCIidcfXzyxl5NIHgttbb/oT0S0b4GJ7j0WE80Nk0aNqdyzc4mM0IiJyNCpM7s2s/+EfoAOwHdgGtPfqREQavH0ZecTcfB2R5AOwKvk0ej9yrc9R1bzMxkUX1e75VhfViojUF5U9xOoh788YYCDwLYELaXsDXwOn12xoIiL+m3vWvYw5tACAXKJIfuNZCGtY97UvT3azDpAZKB9crYtqRUTqiwr/h3LOneWcOwvYBPR3zg10zg0A+gFraytAERG//Pc3H/OTZfcHt1dedi/Nz+zmY0S1Z90Z13I10xjGXD5NHut3OCIiEqLKZu4P6+qcW3Z4wzn3nZnpWeQi0qB9O3sTA/9yKWGBG4WxotXZ9P3XbUdo1XC4ocN46flAudVOf2MREZHQhfLd8koze9bMhpnZUO9JtStrOjAREb+sXw9X/CyGTQRuCbkzohUdPv3ncbEc57AORXfDZLNW5YiI1Buh/E91LbAc+CXwv8AKr05EpMHZtg1GjIDlu1oyjDTeixhF7r9eI65zG79Dq1Xti66nVXIvIlKPHHFZjnMuB3jE+xERabA2bgwk9mu9q4ryouNp9P5/SDnz+Hsod9u2YAbOQebWA+TlxhAZffx8cyEiUl8d8Z3azDaY2frSP7URnIhIbVk4cz3P9nyU778PbEdEwCuvwBnHYWIPEBUFH0eN5Eeas584ti/a5ndIIiISglAuqB1YrBwDXAIk1Uw4IiK1yzl49zdpnDx1PPeRQTqJvBx1DS+/DBdd5Hd0/moZkUHz3F0AZHyzifanpfgckYiIHEkoT6jNKPaz1Tn3KHB2LcQmIlKjtqzI4j8dJ3P+1LNpRgYAf7Vb+HTWbi6+2Ofg6oB9TYuuqs1eoYX3IiL1wRFn7ks9jTaMwEx+fI1FJCJSw/Zm5JN2wyv0f/23XOi2BOt3RbQkb+YsTj5fX04C5LZsD+mBct5aPaVWRKQ+CGVZzkPFyvnABuDSmglHRKTmbFm+j8W3v0LPOQ8ypnBdiX3LO1xA6ntP0ayLlp4cZqkdYFGgHJaumXsRkfoglOT+Z865EhfQmlnHGopHRKRabdgA778PkQ/cx8QNf+RCDpbYvzu8GTt//zg97poYuD2MBMWcVHQ/zNidmrkXEakPQrmv2Wsh1h01MxtpZqvNbK2ZTSlnf7SZzfD2f21mqcX23eHVrzazEUfq08w6en2s8fqMOtIYIlJ/5Oc5Ni/OYP7fFvLWz2Zz5ZXQqROccALccAPM39CMRsUS+z1hTVk89h6a/LCGLndfpsS+HAm9i9bcN83SzL2ISH1Q4cy9mXUFegAJZlb80rImBO6aUyVmFg48CZxLYFXnAjOb7ZxbUeywnwGZzrlOZjYReACYYGbdgYlefG2AD83sJK9NRX0+ADzinHvFzP7h9f33isaoLPbCPMeulRU/j925Yhuxsbi4xiUP2LsXy80Jqb2La4xrFFfygIwMLO9QSO0LmyRCbGyJ/fbjDsjPD619UjOIji6xP2z71lIvsmSj4rsKmreCyMhiHRYSti290kFLtG/TruRTQfPyCP9ha4Wxl2hvRn7bDiX2Wc5BwnZsD619ZCT5rduVbJ+dRfiuHeW22/PtLjbmrw22dzGx5LdsW+KYsL2ZhGVUfO4UH78wLp6CFq1L7AvP+JGwvZmVBF1ULGjSlIJmLUscFrFjK5a1r9LxD8tPakFBYnKJusitG7ED+0OKP69lCoXxCSX2RW38HjuUW3n7QkfhoXz2tDiJA2GNyckh+JPy8UsU7MrEZezG9mQSkbWb6P27ScjeSttDG2hPFu2BA8Qyjn0UFHuLm8mlPMYv2da4C3vGXkOvR6+jX1KTkH4Xx6sWA4tm7lvlbsIVOixMH4LqCucCb+WH386d8/797d+PKygE54rqvPLhA4PvM7FxEBlZdBxgezIpUXG48/JiSEgM3De22HG2OyP019A0qeR7fEFBYPxQ2yc3K1mRl4ft2xta47CwwPjF5eZi2VmhtY+ICLz+4g4ePOJ7ZFBUFC4+8B60d28ku3YB+/djOQcrb+dx0THQuGR+YdlZkFv5e2ywfWwjaNSoZPt9eyEvL7T2cY0hpmQ6aHsyoaAgtPbxTQL33C3efndGhedamfYN/Nw7uGU/Gat3hRxPyeCcK/cHGAO8AGR4fx7+eRw4raJ2of4ApwLvFdu+A7ij1DHvAad65QhgF2Cljz18XEV9em12ARGlx65ojMpiT6Hl4be9I/5M5dYy1f9mYsjtb2Vqmeq5DA25/UT+XaZ6NZ1Dbj+UuWWq99E45PadWV2iqjH7Qm7rwDVmX4mqk1gVctt9NC5TPZS5IbdfxUllqi9jesjt5zK0TPWtTA25/XQuK1M9lVtDbl/fz71T+aJMdS6RIbfvzncOnGvUyLkLLnDu0Ued2/7VRleb5s6dW6vjVbfCgkK3n9jg73T3ut1+h9QgzJ071x3YX+iWfvSje+flve6pp5y75x7nbrzRuQkTnJuV+ks3P/FctzjuNLciqo9bH97JbQ9r7XZZsttNottHY3eAGJdLpLuM6WVO/+/oHvK/k3N4v0x1Bk1Dbt+DZSWqojkYclsHrikZJapSWR9y24NEl6k+jc9Dbr+e1DLV43g15Pafc1qZ6sk8FnL7VxlXpvpP3B5y+8eYXKZ6Gj8Nuf3t/KlM9fucE3L7nzKtTPUyeoTcXudeyerKzj1goXOh59gVztw7594C3jKzU51zXx7bR4dKtQW2FNtOB06u6BjnXL6Z7QWSvfqvSrU9PEVaXp/JwB7nXH45x1c0RomPS2Y2CZgEkELJ2VARqX7RlJ19yiGGKCqeVcomjm3RHdiZ2IGrhq+nyZBcOnXKIirKAbDqIKxK21BjMZeJJzubtLS0WhuvJqREtqNTXuDJXnNfnEPSWW18jqj+yc4KY8fHO3HzvqfFhuV0yFzOofz19GIfj/Isz/OzEsffyGIG8WlIfUdQ8bewInJ8qmxZzm+ccw8Cl5vZZaX3O+duqeLY5X2360I8pqL68q4hqOz4UOPAOfc08DRASlhrt4tmZRqVKzaO5qVW1RzKSmBnbovK23lRRcQ1omXJb83YvyeJHXmtQhs+IYbWpRZRZe5qyfaCyr82PPxLSUiKom3JVTn8uKMN2a5se1fOr7Jli3Byip1ljQqNrT+0K3NcRe3btTEOFPtbbZEXSfqODpXGDoDBAYujY6kbnyTlxLB5Z2jXg2dEpnBiqTwmfn9jNu06sdyzpvSShX0xbelc6q8pem8iGzI7hxI+hxq3okvp0ySjORv2nlT++KUrE1rQrdRpmrujDeuyu4Y0fmyzJHo0LVm/N70Da3O6H7E9QFKrJvQqtepl+8aTsEPR5f+rK6bQIjipUwzhLQLf+h7+WbToSiJjwqFpUyw5icgWTYlunUR8p5a0PKUjjTs04yQzTgKGhBRlzUpLS2PYsGF+h1ElixJSYVcguW/tmnJqPX89tWXrpnwW/ukD7K1ZDPjhbdpS/hN+O7G2TN1Omoc8TpTlExVZtLrADHJyG5FVGFiuUfx9oeR7RKAcGxdO06iitmawLzMRc2XblPceHd8knORi7/HRDjIyk8scV5GmiVZiZURCQRgZe0Nrf4hokkutbIjLiyAjK7T2WWFNSS61qibmUBQZ2aG1PxCRQHKp97iInBgyDoQYf1Q8yd6qmry8PCIjI+FAIzJyQmtfGB1Hcqn8Ij+7MRmHQmsfFhtLcslVu+RkNSEjL7T2UXHRJJfKD7L3JpJREFr7RvGRJEeWrNuTmYRzoVwOepydey705UbgLT8pd4fZT5xz/zGzq8sdx7kXj2qksv2fCtztnBvhbd/h9funYse85x3zpZlFAD8AzYEpxY89fJzXrEyfwJ+BnUArb3Y+OHZFY7iKfjFAly5d3OrVq6vy8qUBagiJnFSvhnBO3D1hJdNnRpBOCg8+HsvkyX5HVHc5B599BktveZYx395DO8q5vqiYbGvMe6k/592zptKqFbRqBcnJkPLjNyQe+pGopnFEJ8URkxxHREIc4Y2iiYiJICI6nMjYCMKjIwiLiii5bljqpYbwXiHVq/g5YWaLnHMDQ21b2bKc/3h/VimJr8QCoLN3W82tBC6QvbzUMbOBq4EvgfHAx845Z2azgX+b2cMELqjtDMwnMBVRpk+vzVyvj1e8Pt+qbIwaes0iIvVKVJ9urJ0ZKG/WDXMqNG8e/N//wfz5cBM53Fwqsd8b3pRNqcNwJ5/CtrbxnH3TWBq3b8k4M8aV6a1/mRoRkVBVtiznP5SzPOUw59yFVRnYm0G/mcAFreHA88655WZ2L4ELB2YDzwH/NLO1wG4CyTrecTOBFQQerHWTc67Ai7tMn96QtwOvmNl9wGKvbyoaQ0REoH3RDXOU3Jdjxw64/XZ4sdg02L+4kr9wGzmR8Wwe+lPa33IRTc8/hd7enT0y09KI7hDa0koRkaNV2UOsptb04M65d4F3S9XdWaycA1xSQdv7gftD6dOrXw8MLqe+wjFERI53Su4r9tmrP7D9qtv4b+5fgECyHh0NE65O5IezP6HjRX1oWupWwiIiNa2yZTmfHC57D3zqSmAmf7VzruKbrIuISIPRvj1EkUs7ttDq+23AmX6H5LvCQvjnjfM4/+mLOIOdxLGb0bzNuHHGQw9Bhw5QzlySiEitqGzmHgAzGwX8A1hHYE17RzP7uXNuTk0HJyIi/mrbMp9sGhNJPuyG3H05RDc5fmejCwrgseGzufGTCcQSeBjhKN7lq4fmcfKtdeEeTSJyvAvlEvuHgLOcc8Occ0OBs4BHajYsERGpCyJjI/gxvOhJyTsWVX4HmIYsPx+eP+1ZfvnJ2GBinxnZnIxXPlBiLyJ1RijJ/Y/OueI34l0P/FhD8YiISB2zq1HRwvvdizf5GIl/8vPh70Nf4fr51xNOIQA/NjmRxt/OI3nCOT5HJyJS5IjLcoDlZvYuMJPAmvtLgAVmdjGAc+6NGoxPRER8lpXUAbK+AGD/yuPzqtq/XZrGz+cVPfZlU/MBtF/2LtbyCA8kFBGpZaEk9zHADmCot70TSAJ+QiDZV3IvItKA5bVqD96Eff6642/mfsady/npmxcRTeBeEj8kdqX9yvex0o+oFBGpA46Y3Dvnrq2NQEREpG6y1A7wdaAcvvX4mrn/5rP99P3DOBLZC0BmTCtaLJqjxF5E6qxQ7pbTEZgMpBY/vqoPsRIRkfohtmuHYDlu1/Ezc5+VBbdftplnOQjAgbA4Yj56l7ATUv0NTESkEqEsy5lF4Cmu/wHvKiIRETluJPYuuqA2Kev4Se5/9Sv4cGs3+vAtz0bcyJn3j6DFaf38DktEpFKhJPc5zrnHazwSERGpk1oOKkruW+ZtwRUUYuGh3Gyt/vrgA3jmmUB5L4kcmvZvWlzub0wiIqEI5d35MTO7y8xONbP+h39qPDIREakTElLiyaQpADHkkrm6Yd8NOS8PbrmlaHv8eLjscgMz/4ISEQlRKDP3vYCrgLMpWpbjvG0REWngzOCHmA7k54SziQ40WpNFUvdWfodVY2ZO+YbkVQeA04mPh7/+VXm9iNQfoST3Y4ETnHOHajoYERGpm34zdD5vvxcJwJsOuvscT0358YdCuj36cz5nIS8zkb3/+xCtWrXxOywRkZCFsiznWyCxpgMREZG6q21qZLC8uQHfDfOdS6fRv3AhAGNtFv9zpea1RKR+CWXmviWwyswWALlenXPOjam5sEREpC5pX3RNbYNN7lcuyOb8z34b3E6feBudTkr1LyARkWMQSnJ/V7GyAacDl9VMOCIiUhd1KLrVfYNN7hff8BSXswOAXdFt6fTM7T5HJCJy9EJ5Qu0nZtYXuBy4FNgA/KOmAxMRkbqjQ8schvIVHdhEtwU5wM/9DqlarV9+kLO+mRrc3vd/d9IsLs7HiEREjk2Fyb2ZnQRMJDBLnwHMAMw5d1YtxSYiInVEh4Q9pBF4+8/c1JSGltx/df1zXM4PAOyMbssJd1/tc0QiIsemsgtqVwHDgZ845053zv0VKKidsEREpC5p1bsFuUQB0NRlkrsry+eIqs/mtYc448sHgtt7J/0GoqN9jEhE5NhVltyPA34A5prZM2Y2nMCaexEROc5ERoexLaLoqtodCxrOwvvPJ71EO9IB2B3Zgk5/vs7niEREjl2Fyb1z7k3n3ASgK5AG/B/Q0sz+bmbn1VJ8IiJSR2TEFSX3md82jOR+317HwLS/BLd3/vTX0KiRjxGJiFTNEe9z75zb75yb7pwbDaQAS4ApNR6ZiIjUKdlJRbfMObByk4+RVJ/X/r6TXS4ZgKywJpz08A0+RyQiUjWhPMQqyDm32zn3lHPu7JoKSERE6qb8NkUz9wUb6v/MvXPw8L9aMIR59OMb5v30KaxJvN9hiYhUyVEl9yIicvwK61g0cx+xrf7P3H/+OSxfHiivievHqY9N9DcgEZFqoOReRERCEtu1KLmPy6j/M/d/+1tR+coroUkT/2IREakuSu5FRCQkTfsULctJzq7fM/c7thXwxmuFwe0bb/QxGBGRaqTkXkREQtJqULtguWX+Vlxevo/RVM2Xv3mT7/K78H88zHkDd9Onj98RiYhUjwqfUCsiIlJcQoto5oYPZ19BHJtpz2Vbc2iW2tjvsI5J4uyX6MxaHuZXLGm7F7jH75BERKqFknsREQmJGdzS7UO++y6wfVoGNEv1NaRjsjxtJ0Oy5gS3O917tY/RiIhULy3LERGRkLUvWnbP5np6Te2a+2cSSWBJ0epmp9G49wk+RyQiUn2U3IuISMjqe3JfWAjtPv1XcPvQJVf6GI2ISPVTci8iIiGr78n9gpfXMuDQVwDkEUHXOy/1OSIRkerlS3JvZklm9oGZrfH+bFrBcVd7x6wxs6uL1Q8ws2VmttbMHjczq6xfC3jcO36pmfUv1td/zWyPmb1d069bRKS+6xKXzh38kb9zA6e++//8Dueo7Xh4erC8MvV8Ilsl+xiNiEj182vmfgrwkXOuM/CRt12CmSUBdwEnA4OBu4p9CPg7MAno7P2MPEK/5xc7dpLX/rC/AFdV2ysTEWnAUhvv4o/8jht4igEbX/M7nKOSm+PosaRoSU70/2hJjog0PH4l92OAF73yi8BF5RwzAvjAObfbOZcJfACMNLPWQBPn3JfOOQe8VKx9Rf2OAV5yAV8BiV4/OOc+ArKq9+WJiDRMzQcWPaW2Zc5mcM7HaI7O/Ke+4cTCtQBkWTwn/eonPkckIlL9/EruWzrntgN4f7Yo55i2wJZi2+leXVuvXLq+sn4r6ktERI5C626JZBG4t30jDpCzNcPniEK374U3guXvu43BGsX6GI2ISM2osfvcm9mHQKtydv0u1C7KqXOV1B9LXyEzs0kElvTQvHlz0tLSjqa5HAeys7N1XkgJDfWcaBXenq4FKwD48Pl3aHxmhyO08F9BASxaHUs/WtOG7fxwWl/f/m4a6nkhx07nhJRWlXOixpJ759w5Fe0zsx1m1to5t91bHvNjOYelA8OKbacAaV59Sqn6bV65on7T/397dx5lVXUnevz7Y5RWQNESkUENIgoOmKBiQBsiDolxymDMSxTfitoajW3S/d6zMzwzvF5Jd7oz2Cuxg8QktrFj4ojzgJYao0aM4IQGRGQWBWQQEaH2++Meqq5lDbeGe09x653u1R8AABf7SURBVPtZqxZ7n7vPPr+rv3Xvr07tuy8wvJlzSpJSmg5MBxg9enSaPHlyyyeo26mtrcW8ULFqzYknBu4DawrF/dCtAzh8B3iOjzwCV2yezLf5Op8Y9CQzfzKOHjvnc+e+WvNC7WdOqLGO5ERey3JmAtt3v5kG3NbEmHuBEyJit+yDtCcA92bLbTZExIRsl5xzis5vbt6ZwDnZrjkTgHXbl+9Iktpm0x4Nd+o3vbRj7Id5yy2FfxM9GPbZo3Mr7CWp3PIq7n8AHB8R84Hjsz4RMT4iZgCklNYA3wOeyn6+mx0DuAiYASwAXgHubmle4C5gYTb+auDL2wOJiEeBPwDHRcTSiDixLM9YkqpE3fCG4r7u1ddyjKQ0KTUU9wBnnJFfLJJUbmVbltOSlNJq4Lgmjs8GzivqXwNc08y4g9swbwIubiaWY9oSuyR1d31HjShsNgz0WdH1i/s5c+C1LMyBA2HKlHzjkaRyyqW4lyTtuAYc0nDnvv/arr8sZ9k3r2Imd3Ezn6LP1FPp08cvrpJUvSzuJUltsuf4EfXtmne6/p37vR7+HeN5hFO4g6f3mA6cn3dIklQ2FveSpDYZfPje/KLHRbxaN4LX0j5MX5/oP6CpHYfzt+z5tYx7+zEA6ggO/N+n5hyRJJWXxb0kqU169O7Jv4/8OfPnF/r/tAgOPTTXkJr18n/cx1C2FdoDjuCgDw3OOSJJKq+8dsuRJO3A9tuvob1oUW5htKrH3XfWt986+hM5RiJJlWFxL0lqs+Li/tVX84ujJVs213Hwkrvr+0POOznHaCSpMizuJUlttiMU989d8xR78CYAb/QczD5nfDjniCSp/CzuJUlt9uG62dzKaTzHwZx607TWT8jBmuvuqm8vGPVxoqdveZKqnx+olSS12dA93mUMMwGY92afnKNp2pBnGor7Pqe73l5S9+BtDElSm+159Mj69t6bXyHVpRyj+aCls1dy8ObZAGylJwddekLOEUlSZVjcS5LabPcxg9nIzgAMZD1r5q/OOaL3m3vtXDbTF4AXd5vE3wwZmHNEklQZFveSpDaLHsHyvh+q77/+p1dyjOaDfrX8RHZnNacwkwVnfj3vcCSpYizuJUntsnq3hqU56+d0neK+rg4eegg2sTN3cAoHXOKSHEndh8W9JKld3hnSUNy/N6/rFPdz5sCaNYX24MEwdmy+8UhSJVncS5LaZ2RDcd/rta5T3M+a1dD+2McgIr9YJKnSLO4lSe2y09iG4r7/qq5T3Pf+9dWcyD38DW9z3HF5RyNJleU+95KkdtnjqIbivmZD1yjut2zcwvkvXsZlbGILvXn90MXAXnmHJUkV4517SVK7DJ84gq30BGDwthVseWtTzhHBvF89wc4U4ljRazjDj7Cwl9S9WNxLktql34DeXLr79RzLwwxlKa+u7Jd3SKz9wwP17df2n5pjJJKUD4t7SVK7/XXcmTzKsSxnKH+dn/8nV3ef01Dc9zzR4l5S92NxL0lqt1GjGtrz5+cXB8DGZes4aMOfAagjOODvpuQbkCTlwOJektRuXam4f3n6w/RiGwAv7XQ4NQftkW9AkpQDi3tJUrttL+77s56Nc/PdMWfT7Q0b3K8Y65IcSd2TW2FKktptTN9XWMUEaniTV5/aH8jv9v3e8xrW2+98msW9pO7JO/eSpHYbfsRe1PAmAMO2LmLzxq25xLHm+eWM3PwiAJvpy5gLJuUShyTlzeJektRufXbbmZU99wagN1tZ8uiiXOKY/4uGJTnPD5jIgMH5b8spSXmwuJckdciKgQfWt1f/cV4uMdyydgoX8XNu4lOsPPr0XGKQpK7A4l6S1CEb9m4o7jfPeSmXGG56chj/yUV8hpvY5fKv5BKDJHUFFveSpA6pO6ChuO+9oPJ37hcvhgULCu1+/eDooysegiR1GRb3kqQO6ffhg+rbA1dW/s79rIbl9kyaBH37VjwESeoyLO4lSR2y57ENd+6HbZgHKVX0+k/fvhwoXPO44yp6aUnqcizuJUkdMuyooWxgFwB2TW+xadGqil07vbeV7996IEsYzq84l+MnbqrYtSWpK7K4lyR1SO8+wat9G+7er3iocktzFv3hKfqnDQxjGVPjQQ6b4BaYkrq3XIr7iBgUEfdHxPzs392aGTctGzM/IqYVHf9IRDwXEQsi4sqIiJbmjYIrs/HPRsSHs+PjIuLxiHghO/65Sjx/Sao2b+5+IHUEC9mPlQs2Vuy6K69vWHD/8vCp9OwVFbu2JHVFed25vxyYlVIaBczK+u8TEYOAK4CjgCOBK4p+CbgKuAAYlf2c1Mq8Hy8ae0F2PsAm4JyU0thsjp9ExK6d+DwlqVuYdfKP2Jm3GclCHt7l5Ipdt/8TD9S3t02ZWrHrSlJXlVdxfxrwm6z9G6Cpbxw5Ebg/pbQmpbQWuB84KSKGAANSSo+nlBJwbdH5zc17GnBtKngC2DUihqSU/ppSmg+QUloOrAJqOvWZSlI3MOIjNWymsCTmxRcrc81t69/mgNV/qu/v+yU/TStJeRX3g1NKKwCyf/dsYsxQYElRf2l2bGjWbny8pXmbm6teRBwJ9AFeacfzkaRubcyYhnalivsFv36UPrxXuGavQxg1aXBlLixJXVivck0cEQ8AezXx0DdKnaKJY6mF4+2Zq/Bg4a8B/wVMSynVNTlBxAUUlvRQU1NDbW1tK5dUd7Nx40bzQu/TnXJi3bpewCQAXnxxGw8++Cg9ynz7aMuMGxmdteftfTSrHq4t7wU7SXfKC5XGnFBjHcmJshX3KaVmFz9GxOvZspgVWWHd1L5pS4HJRf1hQG12fFij48uzdnPzLgWGN3VORAwA7gS+mS3Zae75TAemA4wePTpNnjy5uaHqpmprazEvVKy75cQhg5YxZM3zHPTuPEb1Pp3hx+xb1uu98upl9e0Bnzpth/lv3d3yQq0zJ9RYR3Iir2U5M4Htu99MA25rYsy9wAkRsVv2QdoTgHuz5TYbImJCtkvOOUXnNzfvTOCcbNecCcC67BeAPsAtFNbj/6GTn6MkdSs/i0u4l5P4CV/lzdseK+u1Ni9exciNcwF4j14c9HfHlvV6krSjyKu4/wFwfETMB47P+kTE+IiYAZBSWgN8D3gq+/ludgzgImAGsIDCGvm7W5oXuAtYmI2/GvhydvxM4Fjg3IiYk/2MK89TlqTqtmlEw173m+eUd6/7V65+sL49Z6ejGXbgLmW9niTtKMq2LKclKaXVwAe2NUgpzQbOK+pfA1zTzLiD2zBvAi5u4vh1wHVtDF+S1IQeYw+CZwrt3gvmlfVac+f1YRPj+QhP8/rB7pIjSdv5DbWSpE6x60cbtsypef25sl7rp0s+xZE8xZ6sou7CL7d+giR1Exb3kqROse/JY9mWva0M3zyfbevfLst11q6F2bML7TWxB5PO8OtJJGk7i3tJUqeoGdGPhT0PAKAHiWX3vVCW6zz0ENRlmxaPHw+DBpXlMpK0Q7K4lyR1mmU1h9W333hgblmucf/9De3jjy/LJSRph2VxL0nqNJtGNRT3W58uQ3H/7rtM+/UUvs0VfJTHmHpca99hKEndSy675UiSqlPfIw+DRwvt/gs7v7hfcdOfmLC5lgnU8sX4LcMmLuj0a0jSjsw795KkTjPkxIY798PWPgupc++sr7yuYU3OvGEn0Ldvp04vSTs879xLkjrN/sfuzZ85kgWMZHYazxVvvsfAmj6dNv+AJ+6rb2877oROm1eSqoXFvSSp0/TpG1xw2JPMzVbknPI8TJnSOXPXrXqT/db+BYCt9GTkeZ00sSRVEZflSJI61fjxDe3t+9F3hkVX3U0PCst8nu49gbEfHdh5k0tSlbC4lyR1qnIV9+/8/vb69mtjTyai8+aWpGphcS9J6lRlKe63bGGfl+6p7w74wimdNLEkVReLe0lSpzrkEDi/5y/5b87ivoUjWfenjn9T7epbHmGXug0AvMq+TPjS2A7PKUnVyOJektSp+vaFL+wyk7O4gZEsZOnNf+7wnCtm3FHffmboKey6m2tyJKkpFveSpE63ZvSE+vaWB//Y4fnimWfq23Wf+GSH55OkauVWmJKkTtfv+GMgu2Ff8/KjHZrr3XfhqHdq2Z+5nMydnHPJ33ZChJJUnbxzL0nqdKO/eASbKXx97LBN89m2bGW753rkEXh7UzCXcdww8hsccIhfSytJzbG4lyR1un1H92VOnyPr+4v/+7F2z3VHw3J7PvlJ3AJTklpgcS9J6nQRsOJDk+r76+5s39KclOD2hu3tOfnkjkYmSdXN4l6SVBY9pxxb3x70l1ntmmPhv97IJ1+9kiEsZ8AAOPbY1s+RpO7M4l6SVBb7nXMM79IHgBHrn2fb4mVtnqPXlT/iSv6epQzjO4ffSl+X20tSiyzuJUllMfbInXmiT8Ot9sXT72lh9AelRa+xz/LHAaijBwedP6mVMyRJFveSpLLo0QOWHXISAO/RiyWPLW7T+cu+f219u7bnVCZ/Zo9OjU+SqpHFvSSpbPqd/VlO5xZ2ZzXf2Pqd0k+sq6Pv735d350/8VyX5EhSCSzuJUllM/HzI7iN09nAAB5/HNasKe28LbMepWb9QgDWsitjvn56GaOUpOphcS9JKps994Sjjiq0t22Dm28u7bzXv31VffuOXT7PMcfvVIboJKn6WNxLksrqc59raN9wQwknLF7M3o/fWN9df9YF9PDdSpJK4sulJKmszjwTerKNyTzEmQ9cwFu/uqXF8auu+Bk90zYAapnMKd8aV4kwJakqWNxLkspq6FC4cr8f8xAf43yuZvP//efCV8825Y03GPDbn9d3n5z4VUaMqFCgklQFLO4lSWU34CvTeIfCuvm9lj7NtgcfbnLcktkr+et7HwLgBcYw6fufrFiMklQNLO4lSWX36QtruGGnc+v7ay7+JtTVfWDcP/32EA7nL5zNtVx98E/56CTfpiSpLXzVlCSVXb9+sP68r7GVngDUvPwY71151fvGPPwwXH891NGT6zibz/7nVCLyiFaSdlwW95KkijjrW6P4cZ/L6/vxj1+DGTNg3jzWLN7I2Wc3LMU/9VSYODGnQCVpB5ZLcR8RgyLi/oiYn/27WzPjpmVj5kfEtKLjH4mI5yJiQURcGVG4t9PcvFFwZTb+2Yj4cHZ8n4h4OiLmRMQLEXFhJZ6/JHVHe+4JA374LV5gDAC9tm2B88+HMWOYO/Z/sGxJYYecQYPg5z9vaSZJUnPyunN/OTArpTQKmJX13yciBgFXAEcBRwJXFP0ScBVwATAq+zmplXk/XjT2gux8gBXAR1NK47LrXB4Re3fi85QkFTn/kr5896i7eInR7zs+ZePt/ITLAPjlLws77EiS2i6v4v404DdZ+zdAU98rfiJwf0ppTUppLXA/cFJEDAEGpJQeTykl4Nqi85ub9zTg2lTwBLBrRAxJKW1JKb2bjemLy5Qkqax69IBf3LMP/3DUY/wHl7CY4WyiH3dzEnf3PIXrroPTm3pHkCSVpFdO1x2cUloBkFJaERF7NjFmKLCkqL80OzY0azc+3tK8zc21IiKGA3cC+wP/K6W0vKmAI+ICCnf9qampoba2tsSnqu5i48aN5oXex5xo3qXf6cHNN1/GVY/+M2+91ZtDD13HGWcsZejQWqr9P5l5ocbMCTXWkZwoW3EfEQ8AezXx0DdKnaKJY6mF4+2Zi5TSEuDQbDnOrRFxY0rp9Q8MTmk6MB1g9OjRafLkya1cUt1NbW0t5oWKmRMtO/HE4l4/mn7LqD7mhRozJ9RYR3KibMV9Smlqc49FxOvZspgV2TKbVU0MWwpMLuoPA2qz48MaHd9+t725eZcCw5s5Z3u8yyPiBeAY4MZWnp4kSZLU5eS1xnwmsH33m2nAbU2MuRc4ISJ2yz5IewJwb7bsZkNETMh2yTmn6Pzm5p0JnJPtmjMBWJf9AjAsIvoBZNeYCLzcqc9UkiRJqpC81tz/APh9RHwJWAx8FiAixgMXppTOSymtiYjvAU9l53w3pbQma18E/JrC33Hvzn6anRe4C/gEsADYBPzP7PhBwL9HxPblPv+WUnquDM9XkiRJKrtcivuU0mrguCaOzwbOK+pfA1zTzLiD2zBvAi5u4vj9wKFtDF+SJEnqktz6UZIkSaoSFveSJElSlbC4lyRJkqqExb0kSZJUJSzuJUmSpCphcS9JkiRViSjsEqm2iIgN+GVX+qA9gDfzDkJdijmhppgXasycUGPFObFPSqmm1BPz+hKrHd3LKaXxeQehriUiZpsXKmZOqCnmhRozJ9RYR3LCZTmSJElSlbC4lyRJkqqExX37TM87AHVJ5oUaMyfUFPNCjZkTaqzdOeEHaiVJkqQq4Z17SZIkqUpY3LcgIk6KiJcjYkFEXN7E430j4obs8ScjYt/KR6lKKyEvjo2Iv0TE1oj4TB4xqrJKyImvRcSLEfFsRMyKiH3yiFOVU0JOXBgRz0XEnIj4Y0SMySNOVVZreVE07jMRkSLCHXSqXAmvFedGxBvZa8WciDivtTkt7psRET2BnwEfB8YAn2/ixfdLwNqU0v7Aj4F/qWyUqrQS82IxcC5wfWWjUx5KzIlngPEppUOBG4F/rWyUqqQSc+L6lNIhKaVxFPLhRxUOUxVWYl4QEf2BS4EnKxuhKq3UnABuSCmNy35mtDavxX3zjgQWpJQWppS2AL8DTms05jTgN1n7RuC4iIgKxqjKazUvUkqLUkrPAnV5BKiKKyUnHkopbcq6TwDDKhyjKquUnFhf1N0Z8ANw1a+UugLgexR+4dtcyeCUi1Jzok0s7ps3FFhS1F+aHWtyTEppK7AO2L0i0SkvpeSFupe25sSXgLvLGpHyVlJORMTFEfEKhULu0grFpvy0mhcRcTgwPKV0RyUDU25Kff/4dLas88aIGN7apBb3zWvqDnzjOyuljFF18f+5Gis5JyLii8B44IdljUh5KyknUko/SymNBP4P8M2yR6W8tZgXEdGDwhLff6hYRMpbKa8VtwP7Zss6H6BhxUizLO6btxQo/u1oGLC8uTER0QsYCKypSHTKSyl5oe6lpJyIiKnAN4BTU0rvVig25aOtrxO/A04va0TqClrLi/7AwUBtRCwCJgAz/VBtVWv1tSKltLroPeNq4COtTWpx37yngFERsV9E9AHOAmY2GjMTmJa1PwM8mPzigGpXSl6oe2k1J7I/tf+CQmG/KocYVVml5MSoou7JwPwKxqd8tJgXKaV1KaU9Ukr7ppT2pfD5nFNTSrPzCVcVUMprxZCi7qnAvNYm7dWpIVaRlNLWiLgEuBfoCVyTUnohIr4LzE4pzQR+CfxXRCygcMf+rPwiViWUkhcRcQRwC7AbcEpEfCelNDbHsFVGJb5W/BDYBfhD9pn7xSmlU3MLWmVVYk5ckv015z1gLQ03ilSlSswLdSMl5sSlEXEqsJVCrXlua/P6DbWSJElSlXBZjiRJklQlLO4lSZKkKmFxL0mSJFUJi3tJkiSpSljcS5IkSVXC4l6SJEmqEhb3kqQmRcTuETEn+1kZEcuK+n8q0zUPj4gZLTxeExH3lOPaklQN/BIrSVKTUkqrgXEAEfFtYGNK6d/KfNmvA/+vhZjeiIgVETExpfRYmWORpB2Od+4lSW0WERuzfydHxMMR8fuI+GtE/CAivhARf46I5yJiZDauJiJuioinsp+JTczZHzg0pTQ36/9t0V8KnskeB7gV+EKFnqok7VAs7iVJHXUY8PfAIcDZwAEppSOBGcBXsjE/BX6cUjoC+HT2WGPjgeeL+v8IXJxSGgccA7yTHZ+d9SVJjbgsR5LUUU+llFYARMQrwH3Z8eeAKVl7KjAmIrafMyAi+qeUNhTNMwR4o6j/GPCjiPgtcHNKaWl2fBWwd+c/DUna8VncS5I66t2idl1Rv46G95kewNEppXdo3jvATts7KaUfRMSdwCeAJyJiakrppWxMS/NIUrflshxJUiXcB1yyvRMR45oYMw/Yv2jMyJTScymlf6GwFOfA7KEDeP/yHUlSxuJeklQJlwLjI+LZiHgRuLDxgOyu/MCiD85eFhHPR8RcCnfq786OTwHurETQkrSjiZRS3jFIkgRARHwV2JBSammv+0eA01JKaysXmSTtGLxzL0nqSq7i/Wv43yciaoAfWdhLUtO8cy9JkiRVCe/cS5IkSVXC4l6SJEmqEhb3kiRJUpWwuJckSZKqhMW9JEmSVCX+PyXOFRE5wknZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 384 ms, sys: 4 ms, total: 388 ms\n",
      "Wall time: 388 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# FD modelling of homogeneous elastic medium\n",
    "# ------------------------------------------\n",
    "dx   = 1.0   # grid point distance in x-direction (m)\n",
    "dt = 0.001   # time step (s)\n",
    "\n",
    "FD_1D_visc_SH_JIT(dt,dx,f0,xsrc,Yl,wl,L,'elast')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we run the viscoelastic modelling run and compare it with the elastic analytical solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nx =  500\n",
      "nt =  502\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvcAAAFNCAYAAACExkA1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8VMXex/HPkAChl9AJvUkJhI7SBSmCooIKIgoqlivKlUeEe6+PclWuiqjIFR8bF1GRYgERsUPAiyIghBI6CBpACD2UEJLM88fZZHdTN8kmG+D7fr32xZw5M3NmwyH8dnbOjLHWIiIiIiIil74ige6AiIiIiIj4h4J7EREREZHLhIJ7EREREZHLhIJ7EREREZHLhIJ7EREREZHLhIJ7EREREZHLhIJ7ERHxYow5Y4ypH+h+iIhIzim4FxG5TBljuhhjfjLGnDLGHDfGrDLGtM+unrW2tLV2b0H0UURE/Cs40B0QERH/M8aUBZYADwELgGJAV+BCIPvlD8aYIGttUqD7ISJSGGnkXkTk8tQYwFo711qbZK09b6391lq7CcAYc48xZpsx5oQx5htjTJ2UisYYa4xp6Epfb4zZaoyJM8YcMMY87srvYYyJMcY8YYw5Yow5ZIy5yVV+p+ubgr97tFncGDPNGHPQ9ZpmjCnucf4JVxsHjTH3penDe8aY/zPGLDXGnAV6GmMGGGM2GGNOG2P+MMZM8mirrqv+KNe5E8aYB40x7Y0xm4wxJ40xr+fvj19EJDAU3IuIXJ52AknGmNnGmP7GmAopJ4wxNwF/B24BKgM/AnMzaWcm8IC1tgzQAljmca4aEALUBJ4C3gHuBNrifEvwlMfc/X8AnYAIoBXQAXjS1Z9+wDigN9AQ6J5BP+4AJgNlgP8CZ4G7gPLAAOAh1/vy1BFoBNwOTHP1oTfQHLjNGJPRdURELmkK7kVELkPW2tNAF8DiBN2xxpjFxpiqwAPA89babdbaROBfQITn6L2Hi0AzY0xZa+0Ja+36NOcmW2svAvOASsBr1to4a200EA20dJUdDjxjrT1irY0F/gmMcJ27DZhlrY221p5znUvrc2vtKmttsrU23lobaa3d7DrehPPhJG2w/qyr7Lc4Hwbmuq5/AOcDTWvffpoiIpcOBfciIpcpV/A+0lobhjPqXgNnBLsO8JprespJ4DhgcEbg0xoMXA/sN8asMMZc7XHumMfc9/OuPw97nD8PlHalawD7Pc7td+WlnPvD45xnOsM8Y0xHY8xyY0ysMeYU8CDOhwtPafuSWd9ERC4bCu5FRK4A1trtwHs4Qf4fOFNtynu8Slhrf8qg3lpr7SCgCrAI5+Hc3DiI86EiRW1XHsAhIMzjXK2M3kKa44+AxUAta2054E2cDygiIlc0BfciIpchY8xVxpj/McaEuY5rAcOA1TiB8N+MMc1d58oZY27NoI1ixpjhxphyrqk3p4HcrlIzF3jSGFPZGFMJZ47+h65zC4BRxpimxpiSrnPZKQMct9bGG2M64MzJFxG54im4FxG5PMXhPFD6i2uFmdXAFuB/rLULgReBecaY0678/pm0MwLY5yr3IM4Ds7nxHLAO2ARsBta78rDWfgVMB5YDu4GfXXWyWrbzL8Azxpg4nA8Duf1GQUTksmKsTftNp4iISOAYY5rifOAo7nrgV0REfKSRexERCThjzM2uaUAVcL5V+EKBvYhIzim4FxGRwuABIBbYgzOv/6HAdkdE5NKkaTkiIiIiIpcJjdyLiIiIiFwmFNyLiIiIiFwmggPdgUtR+fLlbcOGDQPdDSlkzp49S6lSpQLdDSlEdE9IRnRfSFq6JyQtz3vi119/PWqtrexrXQX3uVC1alXWrVsX6G5IIRMZGUmPHj0C3Q0pRHRPSEZ0X0hauickLc97whizPyd1NS1HREREROQyoeBeREREROQyoeBeREREROQyoTn3IiIiUmAuXrxITEwM8fHxge5KoVGuXDm2bdsW6G5IgIWEhBAWFkbRokXz1I6CexERESkwMTExlClThrp162KMCXR3CoW4uDjKlCkT6G5IAFlrOXbsGDExMdSrVy9PbQV0Wo4xpp8xZocxZrcxZmIG54sbY+a7zv9ijKnrce5vrvwdxpi+2bVpjKnnamOXq81irvwHjTGbjTFRxpj/GmOa5e+7FhERuXLFx8cTGhqqwF7EgzGG0NBQv3yjFbDg3hgTBMwA+gPNgGEZBNb3AiestQ2BV4EXXXWbAUOB5kA/4A1jTFA2bb4IvGqtbQSccLUN8JG1NtxaGwFMAV7JlzcsIiIiAArsRTLgr38XgRy57wDsttbutdYmAPOAQWnKDAJmu9KfAL2M884HAfOstRestb8Bu13tZdimq861rjZwtXkTgLX2tMf1SgHWz+9TRERERKRABHLOfU3gD4/jGKBjZmWstYnGmFNAqCt/dZq6NV3pjNoMBU5aaxMzKI8x5mFgHFAM50NAOsaY+4H7ASpXrkxkZKQv71GuIGfOnNF9IV4K5J6wluC4OIpcvEhCxYqgEdFC70r/XVGuXDni4uICdv3rr7+ecePG0bt379S8GTNmsHv3bp544gmeeOIJPvjgg3ztw5w5c1i/fj0vv/wyAElJST79TPbv388vv/zCbbfdBsD69euZO3cuL730kk/Xvf766/nzzz8pUaIEAOPHj+emm26ifPnyNG/enIsXLxIcHMwdd9zBX/7yF4oUydsYcIsWLVixYgWhoaGZlpkzZw7XXnst1atXB2DMmDGMGTOGq666Kk/XBliyZAlbtmxh4kRnlvasWbN4/fXXAShdujTPPfccXbt2zdM1Nm3axGOPPUZcXBxBQUE8/vjjDB48GIB9+/YxatQoTpw4QUREBG+//TbFihXjrbfeolSpUtx5553p2ouPjycyMjJvvyestQF5AbcC73ocjwD+naZMNBDmcbwHJ1CfAdzpkT8TGJxZm0BlnBH9lPxawOYM+nQHMDu7vjdu3NiKpLV8+fJAd0EKmXy/J5Yvt7Z5c2vBeT33XP5eT/ziSv9dsXXr1oBe/80337QjR470yuvYsaNduXJlgfVh1qxZ9uGHH049Pn36tE/1li9fbgcMGJDr63bv3t2uXbs2XX6pUqVS04cPH7a9evWyTz31VK6vk6JOnTo2NjY2V33yh6uvvjr1+l988YVt06ZN6vGvv/5qa9asaWNiYvJ0jR07dtidO3daa609cOCArVatmj1x4oS11tpbb73Vzp0711pr7QMPPGDfeOMNa621Z8+etRERERm2l/Lvw/P3BLDO5iDGDuS0nBhXkJ0iDDiYWRljTDBQDjieRd3M8o8C5V1tZHYtcKbx3JSL9yIiUrBmzICePSE62jkOCYFRowLbJ5FLwJAhQ1iyZAkXLlwAnNHVgwcP0qVLF/bt20eLFi0AiI6OpkOHDkRERNCyZUt27doFwPvvv0/Lli1p1aoVI0aMAJwR9V69etGyZUt69erF77//DkBsbCyDBw+mffv2tG/fnlWrVqXrzxdffEHPnj1p3bo1vXv35vDhwwCsWLGCiIgIIiIiaN26NXFxcUycOJEff/yRiIgIXn31VSIjIxk4cCDgfCM0atQowsPDadmyJZ9++mmufj5VqlTh7bff5vXXX08Z+Ex16NAhunXrRkREBC1atODHH38EYO7cuYSHh9OiRQsmTJiQrk3PnyvA1KlTmTRpEp988gnr1q1j+PDhREREcP78eXr06MG6deuybLd06dL84x//oFWrVnTq1Cn1Z+Zp586dFC9enEqVKgHw4osv8tJLL6Uet2nThlGjRjFjxoxc/ZxSNG7cmEaNGgFQo0YNqlSpQmxsLNZali1bxpAhQwC4++67WbRoEQAlS5akbt26rFmzJk/Xzkwgg/u1QCPXKjbFcB6QXZymzGLgbld6CLDM9QlmMTDUtZpOPaARsCazNl11lrvawNXm5wDGmEYe1xsA7PLz+xQR8a+vvoJHH/XOe+ABqFEjMP0RySVj8u+VmdDQUDp06MDXX38NwLx587j99tvTPcz45ptvMnbsWKKioli3bh1hYWFER0czefJkli1bxsaNG3nttdcAZyrJXXfdxaZNmxg+fDiPuv59jh07lscee4y1a9fy6aefct9996XrT5cuXVi2bBkbNmxg6NChTJkyBXAC4BkzZhAVFcWPP/5IiRIleOGFF+jatStRUVE89thjXu08++yzlCtXjs2bN7Np0yauvTbDWcapgXRERATHjh3LsEz9+vVJTk7myJEjXvkfffQRffv2JSoqio0bNxIREcHBgweZMGECy5YtIyoqirVr16YGsdkZMmQI7dq1Y86cOURFRaVOFwKybPfs2bN06tSJjRs30q1bN9555510ba9atYo2bdqkHkdHR9O2bVuvMu3atWPr1q3p6r700kupPyPP16Npf++msWbNGhISEmjQoAHHjh2jfPnyBAc748phYWEcOHDA69opH478LWBz7q0zh34M8A0QBPzHWhttjHkG5+uHxTjTbT4wxuzGGbEf6qobbYxZAGwFEoGHrbVJABm16brkBGCeMeY5YIOrbYAxxpjewEWcVXRSPkyIiBQ+sbEwfDgkJzvHHTrAwoUK7EVyYNiwYcybN49BgwYxb948/vOf/6Qrc/XVVzN58mRiYmK45ZZbaNSoUepIbMrob8WKFQH4+eef+eyzzwAYMWIETzzxBADff/+9V/B4+vTpdHPrY2JiGDt2LLGxsSQkJKSucd65c2fGjRvH8OHDueWWWwgLC8vyPX3//ffMmzcv9bhChQoZlpszZw7t2rXLsi0g3ag9QPv27bnnnnu4ePEiN910ExERESxbtowePXpQuXJlwPnwsHLlSm66KW8TIdauXZtpu8WKFUv9xqJt27Z899136eofOnQotW5O3iM4zyKMHz8+R/09dOgQI0aMYPbs2RQpUiTDtj0/QFapUoXt27fn6Bq+CugmVtbapcDSNHlPeaTjcebRZ1R3MjDZlzZd+XtxVtNJmz82xx0XEQmUv/0NTpxw0rVqweLFULVq+nIHDsDevZDHh8VELkc33XQT48aNY/369Zw/f95rhDfFHXfcQceOHfnyyy/p27cv7777LtZan5YrTCmTnJzMzz//7DUindYjjzzCQw89xO23305kZCSTJk0CYOLEiQwYMIClS5fSqVMnvv/++yyv6WvffLF3716CgoKoUqWKV363bt1YuXIlX375JSNGjGD8+PGULVs22/aCg4NJThmQAJ/Wcs8s8AYoWrRo6nsNCgoiMTExXZkSJUpw6tSp1ONmzZrx66+/en2jsX79+gw/6Lz00kvMmTMnXX63bt2YPn16uvzTp08zYMAAnnvuOTp16gRApUqVOHnyJImJiQQHBxMTE0MNj0GY+Pj4LO+LvAjoJlYiIpID69fDzJnu4//7v/SB/blz8PTT0LgxDBsGFy8WbB9FcsD9NLj/X1kpXbo0PXr04J577mHYsGEZltm7dy/169fn0Ucf5cYbb2TTpk306tWLBQsWpE5nOX78OADXXHNN6qj5nDlz6NKlCwB9+vRJXZ0FICoqKt11Tp06lbpSzOzZs1Pz9+zZQ3h4OBMmTKBdu3Zs376dMmXKZLqqTtprnUgZBMih2NhYHnzwQcaMGZPuw8L+/fupUqUKo0eP5t5772X9+vV07NiRFStWcPToUZKSkpg7dy7du3f3qle1alWOHDnCsWPHuHDhAkuWLEk9l9l78qXdrDRt2pTdu3enHj/xxBNMmDAh9e8uKiqKhQsX8sADD6SrO378eKKiotK9MgrsExISuPnmm7nrrru49Vb3eLQxhp49e/LJJ84q7LNnz2bQIPeK7zt37vR6DsGfFNyLiFwq/vUvd3rgQBgwIH0Za+Gtt5wg/8AB+OKLguufyCVk2LBhbNy4kaFDh2Z4fv78+bRo0YKIiAi2b9/OXXfdRfPmzfnHP/5B9+7dadWqFePGjQNg+vTpzJo1i5YtW/LBBx+kzsWfPn0669ato2XLljRr1ow333wz3XUmTZrE3XffTdeuXVOn+wBMmzaNFi1a0KpVK0qUKEH//v1p2bIlwcHBtGrVildffdWrnSeffJITJ06k1lm+fLnPP4vz588TERFB8+bN6d27N3369OHpp59OVy4yMjL1Ad9PP/2UsWPHUr16dZ5//nl69uxJq1ataNOmjVcQC85I+1NPPUXHjh0ZOHCg1zKXI0eO5MEHH0x9oDaFL+1mpVu3bmzYsCH1G4Abb7yRe++9l86dO9OwYUO6dOnCokWLsp26k50FCxawcuVK3nvvvdS5+Skf4l588UVeeeUVGjZsyLFjx7j33ntT661atcprOVZ/Mll97SEZa9Kkid2xY0eguyGFTGRkJD169Ah0N6QQ8es9sX07NGvmHpLcuBFatsy47JNPwmTXrMU+feCbb/zTB/GLK/13xbZt22jatGmgu1GoxMXFUaZMmUB347IzduxYbrjhhnRBdGJiIqNGjSI5OZkPP/ywwHdM3rBhA6+88kqG+ymk/Pvw/D1hjPnVWpv9gxIuGrkXEbkUlCoF990HxYo5I/aZBfYAo0e7lwv59lvYs6dg+igiUoj8/e9/59y5c+nyg4OD+eCDD5gzZ06BB/YAR48e5dlnn8239hXci4hcCmrVgrffdh6STfN1fDp16kC/fu7jBQvyt28iIoVQ1apVufHGGwPdjXSuu+466tatm2/tK7gXEbmU1KwJjRplX87zIcGFC/OvPyIiUqgouBcRuRwNHAiuzVNYuxZiYgLbHxERKRAK7kVECrOTJ8GHNaHTqVABevZ0H/u4Y6SIiFzaFNyLiBRmL77orGV/zz3OCjk5cfPN7rTHutIiV7qgoKDUZQsjIiJ44YUXAOjRowfr1q3LcXuLFi3y2on2qaeeynbTqfwwadIkpk6dmmWZadOmeT1kev3113Py5MkcX+u9995jzJgxOa7nad++fT6t9f4vz2WAcfYVkMwpuBcRKayshblz4fRpmDULfv89Z/X793enf/wREhL82z+RS1SJEiW8NieaOHFintpLG9w/88wz+baGeV6lDe6XLl1K+fLlA9ij7KUN7n/66acA9eTSoOBeRKSw+vln2L/fSVeoAH375qx+3boQEeHMv588WbvViuTAQw89RLt27WjevLnXhk4TJ06kWbNmtGzZkscff5yffvqJxYsXM378eCIiItizZw8jR45M3Zl07dq1XHPNNbRq1YoOHTqk2431zJkz3HDDDbRp04bw8HA+//xzwBnVbtq0KaNHj6Z58+b06dMndZOnd955h/bt29OqVSsGDx6cbrnHPXv20KZNm9TjXbt20bZtW6ZPn87Bgwfp2bMnPV3T9urWrcvRo0cBeP/992nZsiWtWrVixIgRAHzxxRd07NiR1q1b07t3bw4fPpzlz23FihWp34i0bt2auLg4rLWMHz+eFi1aEB4ezvz589PVS/tNwMCBA4mMjGTixImpm2wNHz4ccHYYBjJtN2WN+CFDhnDVVVcxfPhwrqR9nYID3QEREcnEZ5+504MHO2vc59T69e4170UEcO/ImuJvf/sbt99+u1eZyZMnU7FiRZKSkujVqxebNm0iLCyMhQsXsn37dowxnDx5kvLly3PjjTcycOBAhgwZ4tVGQkICt99+O/Pnz6d9+/acPn2aEiVKeJUJCQlhzpw51KxZk6NHj9KpU6fU5Rt37drF3Llzeeedd7jtttv49NNPufPOO7nlllsYPXo04OxMO3PmTB555JHUNhs0aEC5cuWIiooiIiKCWbNmMXLkSB555BFeeeUVli9f7rUbLkB0dDSTJ09m1apVVKpUiePHjwPQpUsXVq9ejTGGd999lylTpvDyyy9n+rOdOnUqM2bMoHPnzpw5c4aQkBA+++wzoqKi2LhxI0ePHqV9+/Z069bNp7+rF154gddffz1111dPWbW7YcMGoqOjqVGjBp07d2bVqlV06dLFp2te6jRyLyJSGFkLrhE8AG65JXftKLCXwm7SJOc+9eV1//3p699/v3eZSZOyvWTaaTlpA3uABQsW0KZNG1q3bk10dDRbt26lbNmyhISEcN999/HZZ59RsmTJLK+zY8cOqlevTvv27QEoW7YswcHe46rWWv75z3/SsmVLevfuzYEDB1JHx+vVq5f6IaRt27bs27cPgC1bttC1a1fCw8OZM2cO0dHR6a593333MWvWLJKSkpg/fz533HFHln1dtmwZQ4YMSQ36K1asCEBMTAx9+/YlPDycl156KcNreercuTPjxo1j+vTpnDx5kuDgYP773/8ybNgwgoKCqFq1Kt27d2ft2rVZtuOLrNrt0KEDYWFhFClShIiIiNSf3ZVAwb2ISGG0fTvs3u2kS5XyXvlGRPLVb7/9xtSpU/nhhx/YtGkTAwYMID4+nuDgYNasWcPgwYNZtGgR/Tw3i8uAtTbbHVDnzJnDsWPH+PXXX4mKiqJq1arEu1bIKl68eGq5oKAgEhMTARg5ciSvv/46mzdv5umnn04t72nw4MF89dVXLFmyhLZt2xIaGpqrvj7yyCOMGTOGzZs389Zbb2V4LU8TJ07k3Xff5fz583Tq1Int27f7NCUmODiY5OTk1OPsrpPS58xk9rO7Eii4FxEpjDxH7fv1g5AQ/7R7Bc07Fcmt06dPU6pUKcqVK8fhw4f56quvAGd+/KlTp7j++uuZNm1a6lSRMmXKpJtLD3DVVVdx8ODB1NHkuLi4dEHmqVOnqFSpEkWLFmX58uXsT3nOJgtxcXFUr16dixcvMmfOnAzLhISE0LdvXx566CFGjRqVmp9ZX3v16sWCBQs4duwYQOq0nFOnTlGzZk0AZs+enW3f9uzZQ3h4OBMmTKBdu3Zs376dbt26MX/+fJKSkoiNjWXlypV06NDBq17dunWJiooiOTmZP/74gzVr1qSeK1q0KBczeGbIl3avRJpzLyJSGHkG94MG5a2tlSud1XZ+/tlZUvOJJ/LWnog/TZrk01SaTL39tvPKgbRz7vv165e6HCZAq1ataN26Nc2bN6d+/fp07twZcILqQYMGER8fj7WWV199FYChQ4cyevRopk+fnvogLUCxYsWYP38+jzzyCOfPn6dEiRJ8//33qQ+EAgwfPpzrr7+edu3aERERwVVXXZVt/5999lk6duxInTp1CA8PzzBYT2n7s88+o0+fPql5999/P/3796d69eosX748Nb958+b84x//oHv37gQFBdG6dWvee+89Jk2axK233krNmjXp1KkTv/32W5Z9mzZtGsuXLycoKIhmzZrRv39/ihUrxs8//0yrVq0wxjBlyhSqVavmNVWmc+fO1KtXj/DwcFq0aOH1QPD9999Py5YtadOmjdeHmZtvvjnDdrdv357tz/ByZq6kp4f9pUmTJnbHjh2B7oYUMilP54ukyPU98eefUKOGM8oeFARHjoBr/muuvPsuuB6+Y9AgbWgVYFf674pt27bRtGnTQHejUImLi6NMmTJ+b3fq1KmcOnWKZ5991u9tS/5I+ffh+XvCGPOrtbadr21o5F5EpLD54gv39JkuXfIW2AN4bvjy889O23rQVuSydvPNN7Nnzx6WLVsW6K5IAdOcexGRwqZqVejeHYoUyfuUHICrroKUTWqOHIFsvlYXkUvfwoUL2bRpU7olL+Xyp+BeRKSwufFGiIx0AnGPB+FyrUgR6NjRfbx6dd7bFBGRQknBvYhIYRUa6h5xzyvXOtsAbNgAQEICrFkD33wDv//un8uI+ELP+4mk569/FwruRUSuBB4rTySvW8+//gWVKjkD+v36QZ06zkygDDaBFPGrkJAQjh07pgBfxIO1lmPHjhHih2WP9UCtiMiVwCO4P/Pjev4RaQHvh2pXrnQG+BcsgJtvLuD+yRUjLCyMmJgYYmNjA92VQiM+Pt4vQZ1c2kJCQggLC8tzOwruRUQKi2PH4LrroFcv6NsXevf2X9u1a5NcoSJFThynbNJJ6rKPfdSjZk2oWRN+/RWSkiAxEW67DRYvhv79/Xd5kRRFixalXr16ge5GoRIZGUnr1q0D3Q25TGhajohIYbFsmTMffupUeOopvzadlGyIKuIevW/Hr0ydCvv2wS+/QHQ0NG7snEtMhLvucpbbFxGRS4uCexGRwuL7791pf47aAy+8AN8dcwf3Tw7Zzv/8DwS7vr9t0gSWL3dG8QGOHnXveyUiIpcOBfciIoVFPgX327bBM8/A+9xFf5by/KOHaPXxk+nK1agBs2e7j5csge++81s3RESkACi4FxEpDPbudV4AJUtCp05+adZaeOABZ9nLrTTneIf+PPFKtUzL9+oF993nPp44EZKT/dIVEREpAAruRUQKgx9+cKe7d4dixfzS7MKF8OOPTjo4GN59F4KCsq4zaRKkLNyxfj18/rlfuiIiIgVAwb2ISGGQD1NyLl50Rt5TPPIIhIdnX69mTadsimnT/NIdEREpAAruRUQCLTnZe+TeT8H9hx/Crl1Ounx5eDJlmv25c7BuHWzfnmndRx91P2y7cmXqprYiIlLIKbgXEQm0jRudNe4BqlTxbXg9G8nJ8NJL7uPx46FiReDNN6F0aWe3qunTM60fFgZDhriP//3vPHdJREQKQECDe2NMP2PMDmPMbmPMxAzOFzfGzHed/8UYU9fj3N9c+TuMMX2za9MYU8/Vxi5Xm8Vc+eOMMVuNMZuMMT8YY+rk77sWEUkj7ZQcYzIv66Mvv3RWyQEoUwb+8hfXidq1nadsAbZsybKNRx91pz/+2BnwFxGRwi1gwb0xJgiYAfQHmgHDjDHN0hS7FzhhrW0IvAq86KrbDBgKNAf6AW8YY4KyafNF4FVrbSPghKttgA1AO2ttS+ATYEp+vF8RkUxFRbnTfpqS8+KL7vSDDzrTcgDvbwU2b3YH+hno1MlZ/x7gzBlYtMgvXRMRkXwUyJH7DsBua+1ea20CMA8YlKbMICBl1eVPgF7GGOPKn2etvWCt/Q3Y7WovwzZdda51tYGrzZsArLXLrbUp41GrgbB8eK8iIpn78EPYuRPeeAP69ctzc6tWOS+AokVh7FiPk2FhULaskz55Eg4fzrQdY5ydalO8/36euyYiIvkskMF9TeAPj+MYV16GZay1icApIDSLupkQk5OpAAAgAElEQVTlhwInXW1kdi1wRvO/ysV7ERHJPWOgUSN46CGoXj3PzXnOtb/zTveus6nXShmOB9ixI8u27rzTPUvou+/gyJE8d09ERPJRcACvndGk0rTfD2dWJrP8jD6sZFXefSFj7gTaAd0zKIsx5n7gfoDKlSsTGRmZUTG5gp05c0b3hXgJxD3x558hLF7ckZRfe927ryEy0nuyfNNy5ajqSu9YvJhDWUzNAQgPj2DTpvKuh3R3MGDAoXzo+ZVDvyskLd0TklZe7olABvcxQC2P4zDgYCZlYowxwUA54Hg2dTPKPwqUN8YEu0bvva5ljOkN/APobq29kFFnrbVvA28DNGnSxPbo0cPnNypXhsjISHRfiKdA3BP/+7/uafR9+sDdd3dIX+jHH1Mf4m0CNMmmj6NGwWOPOeno6Ca89FKTLMtL1vS7QtLSPSFp5eWeCOS0nLVAI9cqNsVwHpBdnKbMYuBuV3oIsMxaa135Q12r6dQDGgFrMmvTVWe5qw1cbX4OYIxpDbwF3Git1RfOIlJwTp6EGTOcqTHZjJ77IjERZs50Hz/wQCYFczAtB+Dmm93pH36AU6dy1z8REcl/AQvuXSPoY4BvgG3AAmtttDHmGWPMja5iM4FQY8xuYBww0VU3GlgAbAW+Bh621iZl1qarrQnAOFdboa62AV4CSgMfG2OijDFpP2CIiOSPZctgzBi46iq48cbsy2djyRI45JoxU60a3HBDJgVzGNzXqQNt2jjphAT4Sk8miYgUWoGcloO1dimwNE3eUx7peODWTOpOBib70qYrfy/Oajpp8/2z7pyISE55rm/fsmWem3v7bXf6nnuclXIy1KiRO/3bb07EXqxYlm0PGgTr1zvpr76CoUPz1lcREckf2qFWRCRQ0m5elQf798PXX7uP77svi8IlSzrD8WFh0KOHMz0oG/37u9PffOPsgCsiIoVPQEfuRUSuWPv3w65dTrpECbj66jw1N3Om94O09eplU2HnzmxH6z21bQuVKsHRo87S+FFR7qk6IiJSeGjkXkQkEH74wZ3u2hVCQnLdlLXwwQfu49GjfaiUg8AeoEgR6NvXfez5LYGIiBQeCu5FRALBj1NyfvoJ9u1z0uXLZ/EgbR6lnZojIiKFj4J7EZGClpzs1+D+ww/d6VtvheLF89Rcpnr1cqdXr4Zz5zIvKyIigaHgXkSkoG3ZArGxTrpSJWjVKtdNJSTAggXu4zvv9LFiUpIzNeiNN2DCBJ+qVKsGzZq5r/vTTznrq4iI5D8F9yIiBc1z1L5XL2dCey59/TUcP+6ka9WCLl18rGiMM3/n4YdhyhQ4dsynaj17utPLl+esryIikv8U3IuIFDQ/TsmZM8edHj48B58TihSBxo3dxz5sZgVw7bXu9LJlPl5LREQKjIJ7EZGC9sQTzqtNG++J7Dl0+jQs9thTe/jwHDaQw51qAbp3dwb9AdauhTNncnhNERHJV1rnXkSkoPXo4bzy6LPPID7eSbdqBS1a5LCBXAT3oaHOdTZvdqbtr1njPZovIiKBpZF7EZFLlOcqOTketQfv4H7nTp+rXXONO62HakVEChcF9yIil6CDB91z3o2BYcNy0UijRu703r0+V1NwLyJSeCm4FxEpKBcuQFycX5qaO9fZmRacGT5hYblopH59d3rPHneD2fAM7n/+2Vm2X0RECgcF9yIiBeXrr6FiRejaFWbOzFNTnqvk+Ly2fVqhoVCmjJM+cwaOHvWpWoMGULmykz55ErZvz+X1RUTE7xTci4gUlG+/hcRE+O9/czTHPa2tW2HDBiddvDgMHpzLhozxHr33cWqOMZqaIyJSWCm4FxEpKN9840737ZvrZjxH7W+4AcqVy0OfchHcg4J7EZHCSkthiogUhD17nBdAyZLQuXOumklOTr9xVZ5ce60z/F+/PjRt6nM1BfciIoWTgnsRkYLw7bfudI8eTkCdCz/9BPv3O+kKFaB//zz2a8wY55VDbdtC0aJw8aKzRP7Ro1CpUh77IiIieaZpOSIiBcEzuM/DlBzPte1vvTXXnxHyrEQJZ4PdFKtXB6YfIiLiTcG9iEh+u3gRfvjBfZzL4D4hARYscB/nepUcP9HUHBGRwkfBvYhIfvvlF/f69rVrQ+PGuWrmq6/gxAl3M7mctu83ade7FxGRwFNwLyKS39KukmNMrppJ+yBtEX/9Bn/tNRgxwvm0kIMVczp0cKfXr9dmViIihYGCexGR/OY5375Pn1w1ceoULF7sPs7zKjmeFi1yJvP/9FOO1t+vVcv9EO3p0+7FgEREJHAU3IuI5KfkZGjWzImCixSBXr1y1cxnn8GFC046IgKaN/djH3O51r0x0K6d+3jdOj/2SUREckXBvYhIfipSBGbNgj//dLaVrVAhV8188IE77ddRe8h1cA/Okpgpfv3VT/0REZFcU3AvIlIQgoKgZctcVd2/H5Yvd9JFisCwYX7sF+QpuPccuVdwLyISeAruRUQKOc8HaXv3hpo1/XwBP47c66FaEZHAUnAvIlKIWQuzZ7uP7747Hy6SNri31ueqYWFQpYqTjouD3bv93DcREckRBfciIvmk7nvvwe23OxPmUxaoz6E1a9wL2JQuDTfd5L/+papUyWkcnAj96FGfqxrjPXqvh2pFRAJLwb2ISD6p8sMPzpayd93lbGSVC56j9rfeCiVL+qlznoyBBg3cx3qoVkTkkqXgXkQkP+zcScmYGCddqhT06JHjJi5cgHnz3Mf5MiUnhZ8eqtXIvYhIYCm4FxHJD0uWuNPXXQchITlu4ssv3bN56tSBrl391LeM1K3rTu/bl6OqniP32qlWRCSwggPdARGRy5LndrIDB+aqCc8pOXfd5SyDmW9uucUZva9Tx9klKwdq1oSqVeHwYThzxnlG4Kqr8qmfIiKSpYCO3Btj+hljdhhjdhtjJmZwvrgxZr7r/C/GmLoe5/7myt9hjOmbXZvGmHquNna52izmyu9mjFlvjEk0xgzJ33csIleEw4dh5UonbQwMGJDjJmJjYelS9/GIEX7qW2a6dIExY+CGG6BWrRxVTftQrebdi4gETsCCe2NMEDAD6A80A4YZY5qlKXYvcMJa2xB4FXjRVbcZMBRoDvQD3jDGBGXT5ovAq9baRsAJV9sAvwMjgY/y432KyBVo4UL3cpLdukG1ajluYu5cSEx00tdcA40a+bF/+aBNG3d648bA9UNE5EoXyJH7DsBua+1ea20CMA8YlKbMICDli+lPgF7GGOPKn2etvWCt/Q3Y7WovwzZdda51tYGrzZsArLX7rLWbAM0SFRH/+Phjd3pIzr8QtBZmznQf33WXH/qUz1q1cqcV3IuIBE4gg/uawB8exzGuvAzLWGsTgVNAaBZ1M8sPBU662sjsWiIieRcbC5GR7uNbbslxE2vWwKZNTrpECWep/AKVlJTjp2IV3IuIFA6BfKDWZJCXdlvEzMpklp/Rh5WsyvvMGHM/cD9A5cqVifT8z1sEOHPmjO4LofqSJTRxBcbHmzZl086d7l2ofDRlShOgOgDdux8iKmqHv7uZocZTp1Jh/XqKHznChhkziGvSxOe6yckQEtKV+PggDh+Gzz5bRcWKF/Oxt5cu/a6QtHRPSFp5uScCGdzHAJ5PbYUBBzMpE2OMCQbKAcezqZtR/lGgvDEm2DV6n9G1smStfRt4G6BJkya2Ry7WrJbLW2RkJLovhBUroFgxSEjg+LXX5vieOHXKaSLF009Xp1On6v7tY2ZeeQUOHQKgbWhojtfmj4iA1auddOnSnXOztP8VQb8rJC3dE5JWXu6JQE7LWQs0cq1iUwznAdnFacosBlK2bRkCLLPWWlf+UNdqOvWARsCazNp01VnuagNXm5/n43sTkSvV0087U3PmzOFILn4xz5kD58456ZYtoWNH/3YvS3XquNP79+e4uqbmiIgEXsCCe9cI+hjgG2AbsMBaG22MecYYc6Or2Ewg1BizGxgHTHTVjQYWAFuBr4GHrbVJmbXpamsCMM7VVqirbYwx7Y0xMcCtwFvGmJTyIiK5U7Ys3HEHCZUq5aiatfDWW+7j++93lpksMAruRUQueQHdxMpauxRYmibvKY90PE7QnVHdycBkX9p05e/FWU0nbf5anGk6IiIB9csv3g/S3nlnAXcgD7vUgoJ7EZHCIKCbWImIiNtrr7nTQ4dCuXIF3IE8jtyHh7vT27fDhQt+6JOIiOSIgnsREX+YMgXuvdfZmdbmaDEuAGJivJfHf/RRP/bNV3kM7suUgQYNnHRiImzd6qd+iYiIzxTci4jkVXIyzJgB//kPdO8OS5bkuInXX3eWlwdnkZqICP920SeVKzvzgcBZtufkyRw3oak5IiKBpeBeRCSvVqyA33930hUrQp8+Oap+9iy8/bb7+K9/9WPfcsIYPVQrInKJU3AvIpJXs2e708OGQfHiOar+/vtw4oSTbtAABg70Y99ySsG9iMglTcG9iEhenD4Nn3ziPh45MkfVk5Jg2jT38aOPQlCQf7qWK54r5qR8G5EDaYP7XDx+ICIieaDgXkQkL2bNcubVADRvDm3b5qj6p5/Czp1OumxZGDXKz/3LqXHj4Ndf4ehRePjhHFevU8e9ys/x43DggJ/7JyIiWVJwLyKSW8nJ8O9/u4/HjMnRrlPJyfDss97Vy5TxY/9yo3FjaNMGQkNztYOWMc7Ouik0NUdEpGApuBcRya2vvoI9e5x0+fIwYkSOqi9aBFu2OOlSpeCxx/zcvwDRvHsRkcBRcC8iklueu07dd58TofvIWnjmGffxmDFQqZIf+xZACu5FRAJHwb2ISG5s3QrffeekixTJ8fz0xYvdgW/Jks5U90IjORn+/BPWrIGEhBxXV3AvIhI4Cu5FRHIjKQn693fSgwZ5rzLjQ9Unn3QfP/QQVKni3+7lSfPmUL06dOwIe/fmqnoR1/8uu3bBuXN+7p+IiGRKwb2ISG6Eh8PSpbB9O0yenKOq//mP91z7xx/Ph/7lRbVq7nQu1rovWRIaNXLSycnu9yoiIvlPwb2ISF40aQJNm/pcPC7Oe9R+4kTvWLpQqF3bnc7FWvegqTkiIoGSbXBvjClpjPlfY8w7ruNGxphA7p8oInLJeuEFOHLESYeFFbK59ik8d6n1Q3C/eXMe+yMiIj7zZeR+FnABuNp1HAM8l289EhEpzJYudW9alUP79sHLL7uPn3/emcJS6HiO3OdiWg54r3W/aVMe+yMiIj7zJbhvYK2dAlwEsNaeB3K+s4mIyKVu50644QaoXx9eecWZUO4ja+GBB+DCBee4XTu444586mde+WHkPjzcnd682Xn/IiKS/3wJ7hOMMSUAC2CMaYAzki8icmUZP94J6I8ccUbwi/j+2NJ778G33zppY5yNbXNQvWD5YeS+dm0oW9ZJHz8OBw/6oV8iIpItX/5reRr4GqhljJkD/AA8ka+9EhEpbL7+2lmcHpzo/Pnnfa568KD37rOPPQadOvm5f/5Uq5Y7HRPjrN2ZQ8Zoao6ISCBkG9xba78DbgFGAnOBdtbayPztlohIIXLunPcmVaNGQfv2PlVNmY5z6pRz3KABPPtsPvTRn0qWhMqVnXRiorOhVS6knZojIiL5LzizE8aYNmmyDrn+rG2MqW2tXZ9/3RIRKUSeftq9mVO5cvCvf/lc9eWXYckS9/HMmYX0Idq0ateG2FgnvX8/1KyZ4yY0ci8iUvAyDe6BlDUdQoB2wEacB2lbAr8AXfK3ayIihcCPPzoPz6Z4+WWoWtWnqlFR5Zg40X08dix07+7n/uWXOnUgOtoJ8uPjc9WEgnsRkYKXaXBvre0JYIyZB9xvrd3sOm4BFLb9FEVE/O/ECbjzTveqOL16wT33+FT14EF45pnmqdPVr74apkzJp37mh/ffd75iMLlfHK1FC3d6+3ZISIBixfzQNxERyZQvD9RelRLYA1hrtwAR+dclEZFCIDnZCexTloKsUMFZ8saHYPfkSRgwAE6ccCLZKlXg448vscC2VKk8BfbgrJZTt66TvngRduzIe7dERCRrvgT324wx7xpjehhjurt2qt2W3x0TEQmoNWucFXJSzJzpbCmbjbNnYeBAiIpyjosUgXnzcjVl/bKgqTkiIgXLl+B+FBANjAX+Cmx15YmIXL46dYKFC6FECZgwAW6+Odsq58/D4MGwapU77913oWfPfOxnIecZ3GvFHBGR/JfVA7UAWGvjgVddLxGRK8eNNzoj+E2bZlv06FEYNAh++smd9/DDuxg1qlE+djAfWQvr1zsr5Rw8CGPG5KoZz+UwNXIvIpL/sg3ujTG/4dqd1pO1tn6+9EhEJBCshUOHoEYN73zPp0IzsWsXXH897N7tzps0Cbp3PwBcosE9QNeuztcRACNGOMuA5pCm5YiIFCxfpuW0A9q7Xl2B6cCH+dkpEZECdeGCszFVu3bOjqw+stZZVKZdO3dgbwy8+io89VQ+9bWgGOMsg5ki5cHiHGrYEEJCnPSBA3D8uB/6JiIimfJlh9pjHq8D1tppwLUF0DcRkfwXHQ3dusHs2c7I/aBBzlOx2ThyBG67De6+G06fdvJCQuCTT+Cvf83zQjOFg2dwv39/rpoIDoZmzdzHmncvIpK/sg3ujTFtPF7tjDEPAmUKoG8iIvknPh7++U9o3dqZV5+iVSsoWjTTaufOwXPPQYMGTiCfokEDWLkSbrklH/tc0OrUcadzOXIPmpojIlKQsp1zj3unWoBE4DfgtvzpjohIPjt3Dt55x9lR6uBBd37RovDCC/DYYxkOux8+DG+/DW+8AX/+6X1u9GhnE9vSpfO57wXNDyP3oOBeRKQg+RLc32ut3euZYYypl0/9ERHJH1FRMGsWfPSRs7SNp06dnDUrmzf3yj5zBpYudTagWrzY2WHVU7Nm8PLL0K9fPvc9UPww5x68V8zRtBwRkfzlywO1n/iYl2PGmH7GmB3GmN3GmIkZnC9ujJnvOv+LMaaux7m/ufJ3GGP6ZtemMaaeq41drjaLZXcNEbmM/Oc/MH26d2BfrRq8/jr897/QvDmxsfDdd85KNz17QqVKcPvtzvQbz8C+Rg1nFH/jxss4sAfvaTl+GrnfvNnZ/FdERPJHpiP3xpirgOZAOWOM5yzSskBIXi9sjAkCZgDXATHAWmPMYmvtVo9i9wInrLUNjTFDgReB240xzYChrv7VAL43xjR21cmszReBV62184wxb7ra/r/MrpHX9yci+cRaZ0j9xAln6ZXjx51lWH77zf06ccKJvF2SkuDMoLso9+9/A3C+cm029R3P93XuZfe6EuzrDTt3es/SyUinTvDoo85GVcWK5eebLCT8NHJfpQpUrepMbTp3DvbudVbRERER/8tqWk4TYCBQHrjBIz8OGO2Ha3cAdqdM+THGzAMG4eyAm2IQMMmV/gR43RhjXPnzrLUXgN+MMbtd7ZFRm8aYbTgr/NzhKjPb1e7/ZXYNa226tf1TnP0jnv82uNunN7mpRj9W1xvmlddrxwwaxq5Ov3lABn6uM4yNNa/3yrtxy2RqntqWemwBk0l3v2v8MLsqX4Pn6eHr/4fy8YfSF86giYXNn+RAuWZeeQ+sHkWxpPM+9B5mR7zKyRLVU4+LJp7n4TW+/ewAXm8/m4SgEqnHFc8fYFTU2Oy6DcDFIiFM7+Cs2pry/mud2sLQbU9nUNqma+x4iRq82+p1r1LNjq5k4O70+7kZawlKuMAvxaal5u0v24KPmj7nVa7joUVc+/usTHrsLbpiNxY1+B+vvOt+n0nHPxdl2PW0Vle7ie9q3+suZmHwnim0OLbCp+t/W+tefqrm/XToqO0TqHs6+0nTFvi43gQ2h/bwyv/r5nsJvXDAq98m3RuwBNlE3mj6OvtLed97H0VWp9zFowTbxGz7EFHjCAcuVuHCBeezgLVteZHxfENfImN7kPxhULZtgDOlZMgQJ6BPM2vn8hcW5jx/YK3zySchIdefasLDneAenNF7BfcZOHPG+TkfO+akz571fiUkOJ9UExOdOWFDhnjX//RT56unFGn/X/A8vvZaGDrU+/y778LPP3uXS0mnbevWW2HAAO+8l15yVp/yxejR0Lmzd97f/579J+wUjz+efg+Kv/zF+fToi8mToWZN9/HFi06ffDVjBpQq5T4+cgSeeMK3ukWLOs/9eCgREwMjR/pWv3Jl52ftac0a54EgXzRsCE8+6Z337bfOlEVftG9P8kMPk5jo3I5JSWA+/Zjgr7/EWudWSbaA9b59UtJnrrmOUwOGO+WSnfzyC96m5Maf0tchfTtHut/K0Y4DSE5216+z4CVK74/2ruNR1zNvd4/RHGnUObVucjK0+eTvlDh5MH0o4HnsSm/o9TjHqrfwek89FvyF4IRzXuXSpV2JFX0mE1e2Zmpdk3iR6xeNdvfZVS74zBk2lH7Pt7+TtKy1Wb6Aq7Mrk5sXMAR41+N4BPB6mjJbgDCP4z1AJeB14E6P/Jmu9jJs01Vnt0d+LWBLVtfIqu9hVE35O832NZVx6bI/YqjP9ccxNV32crr7XH8oH6XL3kEjn+t3Z3m67NOU9rl+I3Z4ZZXmtM91LdjSnPbKasx2n+uepnS67O4s97n+dhqnyx7GHJ/rL6d7uuxxTPW5/hyGpcueyjif61/q915XVqTLPksJn+tfxze+Fk19hYRY2769tQ8+aO3cudYeOGDzZPny5XlroDCoUcP9A9q7N9fNjPO4dSdN8mP/LiXJydbu2GG3/O//Wvvll+nPP/207zfrbbelrz9+vO/1H300ff0RI3yv//zz6ev37u17/dmz09dv1sz3+t9+m75++fK+19+yxbvu+fM5+2Vx/Lh3/T17fK6bVKy4XbbM2m++sXbJEmsXLrT23VHv+1z/aLl6duxYa//yF2tHj7Z25EhrX+uywOf6G0pdY8PDnR93kybWNmxo7T8rTvO5/gKGpL8dmOBz/Wk8mv52wPd7bwLPp78d8P3eG8Hs9LcDvt97vfk2/e2A7/deM7Z4ZRUn+3sPWGet7zF2VtNynrDWTgHuMMYMS3veWvto7j5OuC+RQZ71sUxm+Rk9Q5BVeV/7gTHmfuB+gDCqZlBFRPypOBfS5cUTQknOc44SHKciJ6jAcSpymKr8Rr3U104as4+6XnVLlUqkbNmLlCnj/Fm27EWqVr1AtWrnqVYtnurV46lRI56gIPc//507nVdunTlzhsjIyNw3UAg0a9SIolWrcqFKFfb98gvxuZx7X7RoVaApAD/8EEv37j6O8F7iQg4cIHT1aiquXUvZrVspGhdHc+BE69ZsLFnSq2yNkydpnHEz6cQeOkR0mnur3oED1PGxfkxMDLvT1L/qzz+p5mP9vXv38nua+i1PnKCij/W3bdvG4TT12589S6mMi6ezceNGTqRZsrZzYiKZL2Lrbc6cbRyskMiFC0GcP1+EpLOJvOBjXYB27eI4mlSMixeLcPGioXp8AhuzrwY4X8Bcm2a3oGuoz70ZF0/n1Cl47TXvvDOAr0HZ2bPpH2w/5mNduTQYa9PFsc4JY26w1n5hjLk7o/PW2tl5urAxVwOTrLV9Xcd/c7X7vEeZb1xlfjbGBAN/ApWBiZ5lU8q5qqVrE3gBiAWqWWsTPa+d2TVsZj8YoHalBvajm33YftLAqRrNOFa/vVd2lR0/UvrovtQyWTlety2nw7ynJlTf/C0hp49kcL30jcU27szZynW9Ttdc/wXB5+My7nKaJg636EV8Oe8PM7VWf0wRm5R1xwEMHIq4nsSSZd1ZiRcJW7vQq0xWDrS/2WvN8eDzcVTd+E1Gl0rHBgVzsMNN7jIGip2KpdK2lZl3GPfPIDGkNEci+niVKHH0DyrsWZdh7d//+J3atd3/tSaUrcTxZl28ypQ6tJsyv7uDmqw2OooPrcnJhu28ypT5PZqSf+5NVzajds5Wq8+Z2t7zSMrtWU/IiQymZGXQVlytppyvVt/rXPntqyl65oRPGzTF1W/FhdAaXnkVNq8k6IL7a/PUdtI0mFwkmLgGEVwsG+p1qujJWJJKliG5mPdjP0FBzgZSKa/ixb3/LFnS2UypoEVGRtKjR4+Cv3AhtH49tG3rpBs1ytuHpkLvxAln6+KZMzNfHigsDP74wzvv88+dpVgrVYIyZZxpH56v4sWdGzllZ7C002p++cX5QXv+o0n7jzXluEULuOYa73MrVnhvtZy2jmde27beyyABfPONsxGcL7p0ST83a9EiOHmS5GQ4f971iof48862FOc9/oyu2YdDpganT0NcnLORXMc9H3HxbIJTJh4SLsCFhIwv/zmDOEmF1OMiJDGCD3zrOzCXYSRQPPW4NHEM5lOf6iYRxIeM8MqrwmH685VP9c9Qmk/xnpJVh330INKn+keowld4T/dtylY6sCaTGt72U4dIehIc7PzuDQqCtvxKC7uZIkVIfRkDpggUMa6067WneFO2lOyYWqZIEWh7/kdqJ+zxKpfZa2eZtuwvG+5Vv92xbwhNOORb/SpdiC3XMLWuMdD2j0WUTDjpfaunpPH+J7CtVh9Ol67h1WbbHR8RnJyQ7p9KRu1tazSI+BIV3PWTk2gd/YH7Oq5yx44epXLlSgDc/PmoX6217Xz6CyKL4D6/uQLpnUAv4ACwFrjDWhvtUeZhINxa+6DrYddbrLW3GWOaAx/hzLOvAfwANML5kWTYpjHmY+BT636gdpO19o3MrpFV35s0aWJ37Njhzx+HXAYUyElauifc4uOd+DQ52fkPLC7Oe8ryZeHPP52N0d57z3nDGQkN5VjDhoR27erstXBZbGWcsYQEiI11FqhKefb9xAn3s/CZpU+dcuYiXMqKFfMeZEibLlbMGbdKeZ08GUuNGpUpWtT57OZ5Lrcvz+A77Surc768iviy1qLkief/H8aYHAX3WU3L+YIMpqeksNbemIM+ZlQ/0RgzBvgGCAL+4wrCn8GZW7QYZy79B64HZo/jrJCDq9wCnIdvE4GHrR5NhwsAACAASURBVHWGkjNq03XJCcA8Y8xzwAZX22R2DRER8Z+QEGjSBLZtcwK36Gjo0CH7epeU06edByWTPL7ZDAmB3r1h4EDo1QsaNGDzihWX7Ie+M2ecxakOHHAekI6NdZ4lzejPkycD3VtHyrd3Gb1Klcr8XMmSUKKE9zeCmQXsaQP3nAa/kZHRl+w9IYVPVl9UT83vi1trlwJL0+Q95ZGOB27NpO5kYLIvbbry9+JeUcczP9NriIiI/4SHO8E9ODvVXnbBfePGcO+9ziYIrVvDAw/AHXc4U2wuAXFxzkqy+/ZBTIw7iPd8nT5dMH0pUwbKlYOyZZ1XmTI5/zMlcA/ybWEskctGpsG9tTZ1vTzXhk9X4Yzk77DWZjKLTURELitxcc40k99/d5YLnDYt2yqZadkSFixw0pf8TrWJibB9e/rlGJ9+2lk39brrCuWUm9OnnQ9YW7fCnj3OngMpr9hY/16rSBHn8YHKlSE0FCpUcF4VK2acTjkuXz4wz8iIXC6y/edjjBkAvImzRKQB6hljHrDW+vbkh4iIXLqSkpydu8CZo/Dqq7kOWj13qt2U/XYJhde5c87WxatXw5Ytzg5dKWrUcF4BZq3zbOyaNc4zttHRTkCf9hnenCpWzHl7NWtC9erOBmWVK2f8Z8WKmpstEgi+fDZ+Gehprd0NYIxpAHwJPj7WLSIil65y5Zw5DnFxzjIlR4860VsueC6usmmTE4AWwsHtrB07Bjfc4Gz2BM40nC++CPgbOXsWVq6En35yAvq1a50HVHOiaFGoWxfq1YNatZwAPu2rUqWAv1URyYYvwf2RlMDeZS+QwTqMIiJy2TEG6tRxRqjBmZ6Ty+C+Th3354Tjx51VEwvBILfvjh+H7t29d2Ft0SIgn1KsdUbkv/oKvv/eCeovXsy+XtGizqMBzZo5fzZoAPXrO68aNTQ/XeRy4EtwH22MWQoswJlzfyuw1hhzC4C19rN87J+IiARa7dru4H7/fveC9TlkjDM1Z9Uq53jTpksouI+Ph5tucgf2xjjPHzya1/0cfWctREU5zy0sWODMk89KaKjz0HL79s7PvVkzZ2n5NHs/ichlxpfgPgQ4DHR3HccCFYEbcIJ9BfciIpezOh77nv7+e56aCg/3Du779ctTcwUjORlGjoQff3TnzZ4NI0ZkWsWfTpxwLvfWW84zvJkJD3d2Pr36aiegr1dPU2hErkTZBvfW2lEF0RERESmkatd2p/fvz1NTng/VXjIr5jz9NMyf7z6eMqVAAvu9e+HFF53NbjPaE6tsWWf6f79+zlL61arle5dE5BLgy2o59YBHgLqe5fO6iZWIiFwi/Dhyf8mtmBMZCZP/v707j7KqOPc+/n0EBwQRUFQUHEElooIi4lUiioIj5r2aqHHAN0aCQ9SobzTBiNFr1Diba3jjtMTZxERAkYuAdhziAAqCoggYFQRBgSCozM/9ozbZp5s+3WfY+5zu07/PWr2o2qf2c+pore6nq2tXZRypcuGFcMUVqb7lP/8J114Ljz1W/TwsgFat4KST4Ec/ggEDwqFJIiKZclmWM5JwiuuzwPp0uyMiIg1OgjP3mdvCf/ABrF4dtldskJYuhbPPDovdIUyP33VXamtdvv4afve7sNvo6hqnyXTvDuefH87EatUqlbcXkQqRS3K/0t3vTr0nIiLSMCU4c7/11mG7xU8+Cbu7fPAB7L9/USHTs2JF2P9x7tywafuIEaltJ/PMMyF5X7iw+vUjj4Srr4a+fbV+XkRyk0tyf5eZDQNeAFZtuOju76TWKxERaTg6dAhHhq5dG44x/fZb2HLLgsN17x6Se4ApUxpwct+pU3iI9oYbwtOqKWzts3gx/Pzn8MQT1a8ffDDceiscdljibykiFS6X5H5f4CzgSOJlOR7VRUSk0jVrFtaat2xZfRa/QD16wMiRoTxlStiIpsFq3jw8UJuC11+HH/4QPv88vtahA9xyC5x+uk53FZHC5JLc/x9gd3dfXW9LERGpTHfemVioHj3i8pQpiYVtNNzhnnvgF7+ofvDUoEFhvX3btuXrm4g0frnMC7wLtEm7IyIi0jQccEBcnjo1bCPfYHz2Wdh7MqVOrV0Lt9++JxddFCf2bdvC6NHw0ENK7EWkeLkk99sDH5rZODMbHX2NSrtjIiJSmXbcEdq3D+Xly2HOnPL2p5pf/jJMoR9ySOJ/Vli5MizDee65eO3+AQfA22+H/epFRJKQy7KczMWGBhwGnJ5Od0REpNKZhaU5L7wQ6lOmQJcu5e0TAO+8Ex9W9dZbYbechCxfDiecAC+/HF874wy47z5o0SKxtxERqX/m3t3/DiwDjgceAvoB/z/dbomISIMya1Z4yvPQQ+G004oO1yDX3WceVnXyydCnTyJhv/1248T+8svD6h8l9iKStKwz92a2J3AaYZZ+MfAUYO5+RIn6JiIiDcXatfDkk6G8225Fh2twyf3778Pf/hbXr7kmkbCrVsF//mf1xP5nP5vDrbfukUh8EZGa6lqW8yHwCnCiu88GMLNflKRXIiLSsGSeUjtvHqxbV9SBTjWTe/cyH9J0441xeeBA2G+/okOuXx9OlB03Lr52yy3Qs+dcQMm9iKSjrmU5JwNfAC+Z2X1m1o+w5l5ERJqali1hm21Cec2ajY9SzVPnztCqVSgvWgQLFhTZv2LMnl39FKmhQxMJ++tfV/9jwLBhcMUViYQWEckqa3Lv7s+4+6nA3kAV8AtgezMbbmb9S9Q/ERFpKDIPsPr006JCbbJJ9ZNpy7o05+ab460v+/eHXr2KDjliRAi7wSWXpHYWlohINbk8UPuNuz/m7icAHYGpwFWp90xERBqWzKU5n31WdLjM/e7feafocIVZsgQefTSuJzBr/49/wODBcf2EE+C228q87EhEmoy8Drd29yXu/id3PzKtDomISAOVcHLfIB6qfeihsAE9hN82itwh56uvwl72q6Mz3bt1g8cfL+rxBBGRvOSV3IuISBOW4LIcqD5zP3ly0eEKM3duWCMEcP75RU2vu8M558D8+aG+zTbw7LOw1VbFd1NEJFdK7kVEJDcJz9zvs0+8z/vcuWV6qPaOO+CTT+Daa8M+/kWGGjMmrj/8MOy6a1EhRUTypuReRERyk/DMffPmcOCBcX3SpKJDFqZTp/C0a8uWBYd4+224KuNptMsvh+OOS6BvIiJ5UnIvIiK5SXjmHqpvTPPWW4mELLlVq8JynDVrQr1XL/jd78raJRFpwuo6xEpERCS23XYwfHhI8nfeOZGTpyohub/hBnjvvVDecsvwAO1mm5W3TyLSdCm5FxGR3JjBkCGJhsxM7idNCtvNb5L235Tdw372PXvCmWeGxf8Feued6rP0N94Ie+jwWREpIy3LERGRstl1V9h221D+17/CYbGpmzwZJkyAm26C3r3jrTDztGYN/OQnsG5dqPfpAxddlGA/RUQKoOReRETKxqwMS3Meeywun3wybLFFQWH+8Ad4991QbtECHnigBH91EBGph74NiYhI/tauha+/TiRUSZP7tWvhiSfi+plnFhTm88/DBjsbXHstdOlSXNdERJKg5F5ERHI3dmxYS7PFFomtvy9pcj9hAixaFModOsARRxQU5oorYMWKUO7aFS69NKH+iYgUqSzJvZm1M7PxZjYr+rdtlnaDojazzGxQxvUDzWy6mc02s7vNwnYN2eJacHfUfpqZHZAR63/M7F9m9lzan1tEpNHbYouwx/26dYlth3nQQXF5yhRYvTqRsLV79NG4fPrp0KxZ3iEmToQnn4zrf/yjdscRkYajXDP3VwET3b0LMDGqV2Nm7YBhwMFAL2BYxi8Bw4HBQJfo65h64h6b0XZwdP8GtwBnJfbJREQqWeZe9wkcZAXhgdrddw/l1ath2rREwm7sm29g5Mi4XsCSnLVr4ZJL4vqPfwx9+xbfNRGRpJQruT8JGBGVRwA/qKXNAGC8uy9x96XAeOAYM+sAtHb3193dgYcz7s8W9yTgYQ/eANpEcXD3icDyZD+eiEiF6tgx3tt+/vz45KYiZS7NeeONREJubOzYkOBDWEvTvXveIR54AN5/P5RbtYJbbkmwfyIiCShXcr+9uy8AiP7drpY2OwFzM+rzoms7ReWa1+uKmy2WiIjkY/PNYYcdQnn9+vBkaQIOOSQuv/JKIiE39swzcfmUU/I+gOvrr+E3v4nrv/oV7LhjQn0TEUlIaodYmdkEYIdaXhqaa4harnkd1wuJlTMzG0xY0kP79u2pqqrK53ZpAlasWKFxIdVU6pg4oE0bWi9YAMCUUaNYtv/+Rcds0aIV0BOAiRNX8dJLrxd7+G01tmYNh44a9e8fepN33pkVef6/ue++3fjyy10A2G67lRx44FtUVa3Puy+VOi6kcBoTUlMxYyK15N7dj8r2mpktNLMO7r4gWh6zqJZm84C+GfWOQFV0vWON6/Ojcra484BOWe7JibvfC9wLsNdee3lfLbKUGqqqqtC4kEwVOyb23Rc++ACAHu3aJbLovE+fsAPN11/D4sWb06lTXzp3LjpsbNassDvO7Nmw6670PPfcvGbuP/0U/vrXuH777VswYMD3C+pKxY4LKZjGhNRUzJgo17Kc0cCG3W8GAaNqaTMO6G9mbaMHafsD46LlNsvNrHe0S87ZGfdnizsaODvaNac3sGzD8h0REclT5kO1Ce2Y06wZHHZYXH/55UTCxrp0gY8+gunT4U9/yntJzjXXwKpVoXzQQWGjHRGRhqhcyf1NwNFmNgs4OqpjZj3N7H4Ad18CXA9Mir6ui64BnA/cD8wG5gBj64oLPA98HLW/D7hgQ0fM7BXgL0A/M5tnZgNS+cQiIpVil13ickI75gB8P2MiPJV192bQrRv075/Xbe+/D488EtdvvVUn0YpIw5Xaspy6uPtioF8t1ycDP82oPwg8mKVdtzziOnBhlr70yafvIiJNXgoz9xCW5myQ+Mx9Ea6+Gjx6SuvYY6v/EiIi0tBo7kFERPKTmdx/+WViYXv2DGdkAXz8cWIb8RTlzTerb41/ww3l64uISC6U3IuISH66doWpU2HpUpg8ObGwm22WwpaY7nDGGXD33QX9leHXv47Lp54KPXok0CcRkRQpuRcRkfxsvjnsvz+0aZP3g6n1SXxpzsyZ8Pjj4VjZbt3yOnRrwgR48cVQbtYMrr8+gf6IiKRMyb2IiDQYiT9U+/zzcfmoo2DTTXO6zb36rP1PfhI23BERaeiU3IuISIPRuzc0j7Z6eO89+OqrIgNmJvfHHZfzbc88A5MmhfLmm4etMEVEGgMl9yIikj93+OILeO21fx9olYSWLaFXr7g+cWIRwZYvr76259hjc7pt3Tr4zW/i+kUXQceO2duLiDQkSu5FRCR/d9wRTnw97DC4555EQ2duQz9uXBGBJkyI19h37w477ZTTbY8/DjNmhPJWW8FVVxXRBxGRElNyLyIi+dttt7g8Z06ioQdkHCX4wgvxHvN5K2BJzpo1cO21cf2yy2DbbQt8fxGRMlByLyIi+dtjj7iccHLfs2fYiAfCXvcFrfpxLyi5f/DBsMc+QLt2IbkXEWlMlNyLiEj+dt89Ln/ySVionpDmzaFfxlnjY8cWEGTaNJg/P5TbtQtP6tbju++qb3d55ZXQunUB7y0iUkZK7kVEJH+tWsH224fymjUwd26i4TMn2p99toAAmbP2AwaEjerrMXx4fCruDjuEB2lFRBobJfciIlKYzNn7hJfmHH98fD7Wq6/C4sV5BqiqisuZi/izWL4cbrwxrg8dCltumed7iog0AEruRUSkMJnr7jcsVE/I9tvDwQeH8rp1BSzNGTkSxo8Pa2uOPrre5nfdFe+pv8sucN55eb6fiEgDoeReREQKk+JDtQADB8blUaPyvLlFi3Ai7U03wY471tl0yRK49da4PmxYOLhKRKQxUnIvIiKFKWFy//zz8M03ib8FALfcAsuWhfJee8FZZ6XzPiIipaDkXkRECpNycr/PPuEL4NtvC3ywth5z58Kdd8b13/427NYjItJYKbkXEZHC7LFHWL/StWuY8k7BqafG5SefzOGGZcvCbwFff51T/KFDYeXKUD7gAPjhD/Pvo4hIQ6LkXkRECrPddmFKfcYMeOKJVN4iM7kfOxaWLq3nhgkTwnqedu3gZz+rs+nbb8Mjj8T1226DTfRTUUQaOX0bExGRwpilng3vuScceGAor14Njz5azw0TJoR/162D9u2zNnOHyy+P6wMHQt++RXVVRKRBUHIvIiIN2rnnxuX77guJeVYTJ8blo47K2mzkSPj730O5eXP4/e+L66OISEOh5F5ERBq0H/84PlBq+nR4660sDT/7DGbNCuUWLeCQQ2pttnw5XHxxXB8yJLVHBkRESk7JvYiIFO7LL+Fvfwv7SY4YkcpbbL119bX3d9yRpWHmrH2fPlk3q7/6apg3L5Tbtw875IiIVAol9yIiUrgpU+Dkk+GXv4QHHkjtbX7+87j8l79k2Xlzw3p7yLokZ9Ik+MMf4vqdd4Znb0VEKoWSexERKVzKe91v0KMHHH10KK9fH/5QUI179Zn7fv02irF6NQweHK/Z798fTj89nf6KiJSLknsRESnczjtDs2ahPH9+2BozJVdeGZfvvx8+/DDjxfffh4ULQ7ldO+jefaP7r74apk4N5RYtYPjwsOGPiEglUXIvIiKF23RT2G23uL7hgdYUHHkkHHFEKK9bB1dckbFzTuaSnH79Ntqic8yY6rP9N9wAu++eWldFRMpGyb2IiBSna9e4/MEHqb2NWThoasNs+5gx8Nhj0Ys1k/sMM2ZUX35zzDFwySWpdVNEpKyU3IuISHH23jsuV1srk7wePeC88+L6BRfAzJnAaaeFLXW23bbaw7SzZ4e19cuXh3qnTmFTH51EKyKVSt/eRESkOJkz9ykn9xBm7zt3DuXly0MuP7v3mfDkk2HdfbTe5rXXwo6Yn38e2rZsCc8+C9ttl3oXRUTKRsm9iIgUp4Qz9wCtWsETT8QHW82bBwccEE6ZfXf6Jrz8ijFkCBx+OHzxRWjTogWMHg37759690REyqp5uTsgIiKNXGZyP3NmeNp1ww46KenZMyTrxx8Pq1aFGfwrr6y+o84G22wDTz8Nffum2iURkQZBM/ciIlKctm1h++1DeeVK+Oyzkrxtv37w8kvr2LOLZ21z5JHw9ttK7EWk6dDMvYiIFO+UU+C778IsfsuWJXvbXkvH8eGynzDn4H48supU/rxyIJtsAr17hx1y+vXTXvYi0rSUJbk3s3bAU8CuwCfAj9x9aS3tBgFXR9X/cvcR0fUDgYeAFsDzwCXu7tnimpkBdwHHAd8C57j7O2bWHRgOtAbWATe4+1MpfGQRkcr23/9dnvedOBFbtJDOix7nt5ftwG9vG1iefoiINBDlWpZzFTDR3bsAE6N6NVGiPgw4GOgFDDOzttHLw4HBQJfo65h64h6b0XZwdD+ERP9sd98ninGnmbVJ8HOKiEiaMve3z9gCU0SkqSpXcn8SMCIqjwB+UEubAcB4d18SzeqPB44xsw5Aa3d/3d0deDjj/mxxTwIe9uANoI2ZdXD3j9x9FoC7zwcWAe0T/aQiIpKOhQth2rRQ3nTTsO+liEgTV64199u7+wIAd19gZrXtOrwTMDejPi+6tlNUrnm9rrjZYi3YcMHMegGbAXNq67CZDSbM+tO+fXuqqqrq/5TSpKxYsULjQqrRmEjXdhMn8r2o/K+uXZk6eXJZ+5MrjQupSWNCaipmTKSW3JvZBGCHWl4ammuIWq55HdcLiRVeDH8NeAQY5O7rawvg7vcC9wLstdde3ldbL0gNVVVVaFxIpiY3Jq67DqZMCXvdv/JKOC02TY888u9im5NPbjT/rZvcuJB6aUxITcWMidSSe3fPuvjRzBZGy2IWRIn1olqazQP6ZtQ7AlXR9Y41rs+PytnizgM61XaPmbUGxgBXR0t2RESkECNHhuQeQoJ/2GHpvZe71tuLiNSiXGvuRwODovIgYFQtbcYB/c2sbfQgbX9gXLTsZrmZ9Y52wTk74/5scUcDZ1vQG1gW/QKwGfAMYT3+XxL+jCIiTUspT6qdMyfeT79VKzjooHTfT0SkkShXcn8TcLSZzQKOjuqYWU8zux/A3ZcA1wOToq/romsA5wP3A7MJa+TH1hWXsF3mx1H7+4ALous/Ar4PnGNmU6Ov7ul8ZBGRClfK5D5z1r5v3/BArYiIlOeBWndfDPSr5fpk4KcZ9QeBB7O065ZHXAcurOX6o8CjeXZfRERq07VrXJ4xI933Gj8+Lvfb6Nu+iEiTVa6ZexERqTT77BOXp09P971WroRmzUK5f/9030tEpBFRci8iIsnYc0/YfPNQnjcPliypu30xxoyBr76CUaOq/8VARKSJU3IvIiLJaN68+uz9hgOm0tKmDQwcCFbbbsciIk2TknsREUnO/vvH5XffLV8/RESaKCX3IiKSnP32i8tK7kVESk7JvYiIJCftmfuXX4Zf/QpeeglWrUo+vohII6fkXkREktOjB9xzD7zyCrz4YvLx//xnuOkmOPJIGDYs+fgiIo1cWfa5FxGRCtWmDVxwQf3tCpW5v/3RR6f3PiIijZRm7kVEpHH49FP46KNQbtECDj20vP0REWmAlNyLiEjj8Nxzcfnww2GLLcrXFxGRBkrJvYiIpGPxYvjHP5KLl5ncn3BCcnFFRCqI1tyLiEiyvvsOunWDjz+GTTeF5cvjk2sLtWJF9Qd0ldyLiNRKM/ciIpKsFi3i8po1MH168THHj4fVq0N5331hl12KjykiUoGU3IuISPIOOiguv/lm8fEyl+SceGLx8UREKpSSexERSd4hh8TlV18tLtb69TBmTFzXkhwRkayU3IuISPL69InLr7wC7oXHevddWLgwlNu3h169iuubiEgF0wO1IiKSvP32g1atwoOwn38On3wCu+1WWKwePcL+9s8+C+vWQbNmiXZVRKSSKLkXEZHkNW8O//Ef8MILof7qq4Un9wBdusBllyXTNxGRCqZlOSIiko6aS3NERCR1Su5FRCQdmcn9hAnFrbsXEZGcKLkXEZF0HHIItGwZyv/8J8yenX+Me+6BmTOT7ZeISAVTci8iIunYbDPo1y8cODVkCGyS54+cGTPgootg772hd++wJaaIiNRJD9SKiEh6HnsszN6b5X/vQw/F5U6d8v/lQESkCVJyLyIi6WnVqrD71q6FRx6J6+eck0h3REQqnaZBRESk4Rk3Dr74IpR32AEGDChvf0REGgkl9yIiUjorVuTW7o9/jMtnnRX2zRcRkXopuRcRkXS5w4gRcMwxYRZ+2bK628+cCc8/H8pmMHhw+n0UEakQSu5FRCRdZnDnnWGpzTffwFNP1d3+rrvi8oknQufO6fZPRKSCKLkXEZH0nX12XL7jjuzbWn76KTzwQFy/9NJ0+yUiUmGU3IuISPrOPRdatw7lDz+Ml93UdM01sHp1KB98MPTtW5LuiYhUCiX3IiKSvtatq6+dHzoU1qzZuN2OO8blm28ubH98EZEmTMm9iIiUxqWXQosWoTxtWkjea7rxRhg/Hs47Dw4/vLT9ExGpAEruRUSkNHbaCa6/Pq5fc02of/559XZHHQX33lvavomIVIiyJPdm1s7MxpvZrOjftlnaDYrazDKzQRnXDzSz6WY228zuNgt/t80W14K7o/bTzOyA6PouZva2mU01s/fNbEgpPr+ISJN1ySVhLT2ELTKvuSbM1ouISCLKNXN/FTDR3bsAE6N6NWbWDhgGHAz0AoZl/BIwHBgMdIm+jqkn7rEZbQdH9wMsAP7D3btH73OVmWUs+BQRkUQ1bw7PPRcn+AAjR2bfPUdERPJSruT+JGBEVB4B/KCWNgOA8e6+xN2XAuOBY8ysA9Da3V93dwcezrg/W9yTgIc9eANoY2Yd3H21u6+K2myOlimJiKRv223hxRfDkpyePWHTTcMWmCIiUrRynee9vbsvAHD3BWa2XS1tdgLmZtTnRdd2iso1r9cVN1usBWbWCRgDdAb+n7vPr63DZjaYMOtP+/btqaqqyvGjSlOxYsUKjQupRmOiHocdFr4gJPdNJMHXuJCaNCakpmLGRGrJvZlNAHao5aWhuYao5ZrXcb2QWLj7XGC/aDnOSDN72t0XbtTY/V7gXoC99trL+2rvZamhqqoKjQvJpDEhtdG4kJo0JqSmYsZEasm9ux+V7TUzWxgti1kQLbNZVEuzeUDfjHpHoCq63rHG9Q2z7dnizgM6ZblnQ3/nm9n7QB/g6Xo+noiIiIhIg1OuNeajgQ273wwCRtXSZhzQ38zaRg/S9gfGRctulptZ72iXnLMz7s8WdzRwdrRrTm9gWfQLQEczawEQvcehwMxEP6mIiIiISImUa839TcCfzexc4DPghwBm1hMY4u4/dfclZnY9MCm65zp3XxKVzwceAloAY6OvrHGB54HjgNnAt8D/ja53BW4zsw3LfW519+kpfF4RERERkdSVJbl398VAv1quTwZ+mlF/EHgwS7tuecR14MJaro8H9suz+yIiIiIiDZK2fhQRERERqRBK7kVEREREKoSSexERERGRCqHkXkRERESkQii5FxERERGpEBY2kpF8mNlytB++bGxb4Ktyd0IaFI0JqY3GhdSkMSE1ZY6JXdy9fa43lmuf+8Zuprv3LHcnpGExs8kaF5JJY0Jqo3EhNWlMSE3FjAktyxERERERqRBK7kVEREREKoSS+8LcW+4OSIOkcSE1aUxIbTQupCaNCamp4DGhB2pFRERERCqEZu5FRERERCqEkvs6mNkxZjbTzGab2VW1vL65mT0Vvf6mme1a+l5KqeUwLr5vZu+Y2VozO6UcfZTSymFMXGZmM8xsmplNNLNdytFPKZ0cxsQQM5tuZlPN7FUz+145+imlVd+4yGh3ipm5mWkHnQqXw/eKc8zsy+h7xVQz+2l9MZXcZ2FmzYB7jISDgwAABPpJREFUgGOB7wGn1/LN91xgqbt3Bu4Abi5tL6XUchwXnwHnAI+XtndSDjmOiSlAT3ffD3ga+H1peymllOOYeNzd93X37oTxcHuJuyklluO4wMy2Ai4G3ixtD6XUch0TwFPu3j36ur++uErus+sFzHb3j919NfAkcFKNNicBI6Ly00A/M7MS9lFKr95x4e6fuPs0YH05Oigll8uYeMndv42qbwAdS9xHKa1cxsTXGdWWgB6Aq3y55BUA1xN+4VtZys5JWeQ6JvKi5D67nYC5GfV50bVa27j7WmAZsE1Jeiflksu4kKYl3zFxLjA21R5JueU0JszsQjObQ0jkLi5R36R86h0XZtYD6OTuz5WyY1I2uf78ODla1vm0mXWqL6iS++xqm4GvObOSSxupLPp/LjXlPCbM7EygJ3BLqj2ScstpTLj7Pe6+B3AlcHXqvZJyq3NcmNkmhCW+l5esR1JuuXyveBbYNVrWOYF4xUhWSu6zmwdk/nbUEZifrY2ZNQe2BpaUpHdSLrmMC2lachoTZnYUMBQY6O6rStQ3KY98v088Cfwg1R5JQ1DfuNgK6AZUmdknQG9gtB6qrWj1fq9w98UZPzPuAw6sL6iS++wmAV3MbDcz2ww4DRhdo81oYFBUPgV40XVwQKXLZVxI01LvmIj+1P4nQmK/qAx9lNLKZUx0yageD8wqYf+kPOocF+6+zN23dfdd3X1XwvM5A919cnm6KyWQy/eKDhnVgcAH9QVtnmgXK4i7rzWzi4BxQDPgQXd/38yuAya7+2jgAeARM5tNmLE/rXw9llLIZVyY2UHAM0Bb4EQz+62771PGbkuKcvxecQvQCvhL9Mz9Z+4+sGydllTlOCYuiv6aswZYSjxRJBUqx3EhTUiOY+JiMxsIrCXkmufUF1cn1IqIiIiIVAgtyxERERERqRBK7kVEREREKoSSexERERGRCqHkXkRERESkQii5FxERERGpEEruRUREREQqhJJ7ERGplZltY2ZTo68vzOzzjPo/UnrPHmZ2fx2vtzez/0njvUVEKoEOsRIRkVq5+2KgO4CZXQuscPdbU37bXwP/VUefvjSzBWZ2qLu/lnJfREQaHc3ci4hI3sxsRfRvXzP7u5n92cw+MrObzOwMM3vLzKab2R5Ru/Zm9lczmxR9HVpLzK2A/dz93ah+eMZfCqZErwOMBM4o0UcVEWlUlNyLiEix9gcuAfYFzgL2dPdewP3Az6M2dwF3uPtBwMnRazX1BN7LqF8BXOju3YE+wHfR9clRXUREatCyHBERKdYkd18AYGZzgBei69OBI6LyUcD3zGzDPa3NbCt3X54RpwPwZUb9NeB2M3sM+Ju7z4uuLwJ2TP5jiIg0fkruRUSkWKsyyusz6uuJf85sAhzi7t+R3XfAFhsq7n6TmY0BjgPeMLOj3P3DqE1dcUREmiwtyxERkVJ4AbhoQ8XMutfS5gOgc0abPdx9urvfTFiKs3f00p5UX74jIiIRJfciIlIKFwM9zWyamc0AhtRsEM3Kb53x4OylZvaemb1LmKkfG10/AhhTik6LiDQ25u7l7oOIiAgAZvYLYLm717XX/cvASe6+tHQ9ExFpHDRzLyIiDclwqq/hr8bM2gO3K7EXEamdZu5FRERERCqEZu5FRERERCqEknsRERERkQqh5F5EREREpEIouRcRERERqRBK7kVEREREKsT/AiZslyurap/zAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 292 ms, sys: 8 ms, total: 300 ms\n",
      "Wall time: 295 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# FD modelling of homogeneous viscoelastic medium\n",
    "# -----------------------------------------------\n",
    "dx   = 1.0   # grid point distance in x-direction (m)\n",
    "dt = 0.001   # time step (s)\n",
    "\n",
    "FD_1D_visc_SH_JIT(dt,dx,f0,xsrc,Yl,wl,L,'visc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* How to solve the 1D viscoelastic SH problem by a staggered grid finite-difference approach\n",
    "* Optimize the anelastic coefficients $Y_l$ by the global optimization Differential Evolution algorithm to achieve a constant Q-spectrum\n",
    "* Compared the effect of intrinsic damping for $Q=20$ to the elastic analytical solution"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
